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Overview 

The  m^n  goal  of  our  funded  research  was  to  perform  a  systematic  inves¬ 
tigation  and  evaluation  of  the  utility  of  envelope  equations,  in  particular 
the  nonlinear  Schrodinger  equation  (NLS)  and  its  generalizations,  for  mod¬ 
elling  self-focusing  (SF)  collpase  phenomena  in  transparent  dielectric  media. 
The  relevance  of  this  work  for  the  mission  of  the  Air  Force  is  that  envelope 
equations  are  now  being  used  as  models  for  ultrashort  pulse  eye  damage 
studies  and  it  is  a  priori  unclear  whether  these  equations  are  valid.  Specifi¬ 
cally,  for  the  high  peak  powers  present  in  ultrashort  (sub-picosecond)  laser 
pulses  the  peak  power  may  exceed  the  critical  power  for  spatial  self-focusing 
collapse  in  the  vitreous  humor  (1  MW),  resulting  in  potentially  increased 
retinal  damage  due  to  both  increased  loc2d  intensities  at  the  retina  and  also 
the  occurence  of  light-induced  breakdown  (LIB),  It  is  well  known  that  for 
peak  powers  exceeding  the  criticcd  power  for  self-focusing  the  envelope  equa¬ 
tions  predict  that  the  field  will  collapse  to  a  singularity  in  a  finite  distance, 
thus  clearly  violating  the  assumptions  underlying  the  envelope  approxima¬ 
tion.  On  the  other  hand  other  physical  effects  occuring  within  the  vitreous 
humor,  such  as  multi-photon  absorption  (MPA),  serve  to  oppose  the  col¬ 
lapse  and  act  as  perturbations  to  the  envelope  equations.  The  issue  is  then 
whether  the  perturbed  or  generalized  envelope  .^nations  avoid  the  unphysi¬ 
cal  collapse  and  are  therefore  legitimate  model^  for  eye  damage  studies,  and 
our  work  has  shown  that  this  is  indeed  the  case,  i 
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Our  work,  which  was  performed  during  the  years  1993-1996,  was  success- 
fuly  carried  to  completion  and  has  resulted  in  several  journal  publications 
[l]-[9]  and  presentations  at  national  and  international  meetings.  Basically, 
the  research  performed  is  divided  into  three  sections,  namely 

1.  Fundamental  issues  concerning  arrest  of  collapse  in  transparent  media. 

2.  Development  of  reliable  modeb  for  collapse  and  optical  breakdown  in 
transparent  media,  in  particular  water. 

3.  Nonlinear  effects  at  highly  absorbing  interfaces. 

In  the  following  paragraphs  we  provide  overviews  of  the  work  conducted  with 
details  given  in  the  following  Sections.  Section  (1)  describes  our  fundamen¬ 
tal  studies  of  the  arrest  of  collapse  in  transparent  media  displaying  normal 
group- velocity  dispersion  (GVD).  It  had  previously  been  suggested  on  the 
basis  of  numerical  simulations  that  normal  GVD  can  arrest  the  spatial  col¬ 
lapse  associated  with  the  envelope  equations  by  casuing  pulse-splitting.  To 
verify  this  conclusion  we  developed  a  theory  of  critical  collapse  for  pulses  in 
the  presence  of  normal  GVD  [1].  However,  further  studies  showed  that  nor¬ 
mal  GVD  alone  is  not  enough  to  arrest  collapse  [2],  but  rather  the  collapse 
occurs  at  higher  peak  powers  as  the  GVD  is  increased,  that  is,  it  delays  the 
collapse  as  oppossed  to  removing  it  completely. 

Section  (2)  contains  the  bulk  of  the  work  conducted,  and  deals  with 
the  development  of  reliable  models  for  ultrashort  pulse  propagtion  in  wa¬ 
ter  as  a  model  for  propagation  in  the  vitreaous  humor.  In  particular,  we 
have  developed  a  comprehensive  model  for  propagation  in  water  which  ac¬ 
counts  for  self-focusing  collapse,  normal  GVD,  multi-photon  absorption  and 
avalanche  breakdown.  This  work  was  conducted  in  close  collaboration  with 
our  colleagues  at  the  Armstrong  Laboratory,  Brooks  AFB,  and  resulted  in 
the  transfer  of  an  operational  code  for  nonlinear  propagation  in  the  vit¬ 
reous  humor  which  Weis  later  used  for  studies  of  eye  damage.  In  particu¬ 
lar,  we  maintained  regular  contact  with  Drs.  Richard  Albanese  and  Mary 
Potasek  at  the  Mathematical  Products  Division  and  Drs.  Mark  E.  Rogers, 
Pat  Roach,  Randy  Thompson  and  Paul  Kennedy  at  the  Optical  Radiation 
Division.  This  interaction  was  mediated  in  the  form  of  regular  contract 
reviews,  meetings  and  Lab.  visits  which  are  tabulated  in  Sec.  5. 

The  final  Section  (3)  deals  with  nonlinear  effects  at  highly  absorbing  in¬ 
terfaces.  This  problem  is  relevent  to  the  optical  properties  for  light  incident 
from  the  vitreous  humor  onto  the  retina,  in  partilfular  the  highly  absorbing 
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retinal  pigment  epithelium  in  which  damage  can  occur.  At  present  very 
little  is  known  about  the  nonlinear  properties  of  highly  absorbing  interfaces 
but  it  is  clear  that  the  usual  envelope  equations  must  fedl  as  the  absortpion 
per  wavelength  approaches  unity.  We  have  investigated  nonlinear  refection 
properties  in  this  regime  and  have  discovered  a  new  phenomenon,  the  Non¬ 
linear  Optical  Skin  Eflfect  (NOSE)  in  which  a  short  pulse  reflected  from  the 
nonlinear  absorber  suffers  a  Doppler  shift  resulting  from  a  moving  boundary 
induced  in  the  medium  [8]. 

In  the  remaining  Sections  we  now  provide  more  details  of  each  area  of 
research. 

1  Fundamental  Issues 

Our  work  has  served  to  elucidate  the  understanding  of  self-focusing  collapse 
of  ultrashort  pulses.  In  particular,  we  investigated  the  issue  of  whether 
normal  GVD  can  halt  self-focusing  collapse  as  has  been  proposed,  and  also 
whether  collpase  is  a  feature  of  the  full  solution  of  Maxwell’s  equation.  High¬ 
lights  of  this  area  include 

•  An  analytic  theory  of  the  effects  of  normal  GVD  on  collapse  showing 
its  arrest  close  to  the  critical  power  [1]. 

•  Derivation  of  the  threshold  power  for  collapse  in  the  presence  of  GVD 

[2]. 

•  Numerical  study  of  the  conical  emission  and  spectral  broadening  which 
accompany  collapse  and  their  roles  as  signatures  of  pulse-splitting  [3]. 

•  Development  of  a  numerical  scheme  for  soling  the  full  Maxwell’s  equa¬ 
tions  for  nonlinear  optics. 

•  Theoretical  prediction  of  carrier  wave  shocking  for  femtosecond  optical 
pulses  [4]. 

1.1  Arrest  of  crtitical  collapse  by  GVD 

To  address  the  effects  of  normal  GVD  on  SF  collapse  we  generalize  the 
usual  steady-state  SF  equation  to  allow  for  pulsed  operation.  Then  taking 
A{x,y,t]z)  to  be  the  envelope  of  the  singly  polarized  electric  field  E  with 
carrier  wave  exp{ikz  -  wt),  k  =  and  introducing  the  retarded  time 
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i  =  t  —  k'z,  k!  =  is  the  inverse  group  velocity  of  the  wavepacket,  we 

obtain  the  scaled  propagation  equation 

Az  =  +  t|A|*A  -  i'yAtu  (1) 

where  we  now  understand  t  to  represent  the  retarded  time  t—k'z.  In  (1),  the 
transverse  coordinates  have  been  scaled  with  v^,7  is  =  -2^^"-  We 
denote  by  d  the  number  of  transverse  dimensions.  Here  d  =  1  or  2.  We  take 
the  order  of  nonlinearity  s  to  be  such  that  sd  =  4  so  that  here  s  =  4  or  2. 
The  case  where  k"  >  0(jfc"  <  0)  is  called  normal  (anomolous)  dispersion.  For 
pulses  in  the  100  femtosecond  to  1  picosecond  and  beyond  range  with  beam 
widths  of  the  order  of  millimeters,  dispersion  along  the  pulse  propagation 
direction  is  small  when  compared  to  diffraction.  Therefore,  at  least  initially, 
the  term  'yAu  is  small  and  each  individual  cross-section  of  the  pulse  pursues 
its  own  dynamics  independent  of  the  others. 

In  the  absence  of  GVD  each  time  slice  in  the  pulse  is  assumed  to  collapse 
in  the  self-similar  form 

A(t,  C,  t)  =  pJ  x(0  exp[-i^  +  ir]  ,  (2) 

after  capturing  a  critical  number  of  photons.  The  function  x(0  gives  the 
shape  of  each  cross-section,  and  the  new  independent  variables  are  defined 
as 


i  =  9P  ,  (3) 

T  =  j  g^dz'  ,  (4) 

where  a  =  Qr/g,  and  the  operators  transform  as  dz  =  +  ^r)  and 

A  special  combination  of  g  and  its  derivatives  also  arises  in  this 
analysis  and  is  defined  as, 


/3=l(a,-Ha2).  (5) 

This  new  variable  is  proportional  to  the  difference  of  the  local  curvature  and 
twice  the  square  of  the  loced  scale  length  of  g. 

In  the  absence  of  dispersion,  the  equations 


dr 

da 

dr 


otg  , 


4)8  -  a^  , 


dr 


-aw(/3)  , 
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(6) 

(7) 

(8) 


describe  the  evolution  of  g,  a,  and  ^  on  the  collapsing  filament  where  a  = 
dimension  times  the  photon  number  divided  by  the 
width  of  the  ground  state  solitary  wave  of  the  cross  section  at  t  =  0. 

To  proceed  we  extend  the  above  analysis  of  the  self-similar  collapse  which 
has  been  so  successful  in  the  sign  definite  case  to  help  us  understand  what 
the  effect  of  normal  dispersion  is  on  the  pulse  as  it  begins  to  collapse  in 
d-dimensions.  First  we  are  going  to  assume  that  (t)  the  strength  of  the  term 
with  positive  dispersion,  where  to  is  the  initial  pulse  width,  is  small, 
and  (ii)  the  central  cross  section  of  the  pulse,  which  we  take  to  be  t  =  0, 
has  sufficient  photons  so  that  it  collapses  and  reaches  self-similar  form  (2). 

Using  the  self-similar  form  of  A,  as  stated  in  (2),  the  dispersion  term  is 
rewritten  as, 

dtA  =  -dtZod^A  ,  (9) 

d^A  =  -d^zod,A  +  {dtZofdlA.  (10) 

At  the  peak  of  the  envelope  (10)  simplifies.  First,  dtZo\t=o  =  0  for  a  sym¬ 
metric  initial  pulse  centered  at  f  =  0.  Second,  each  of  the  neighboring 
cross-sections  have  longer  collapse  distances  making  5^2o|f=o  >  0.  Thus,  in 
the  neighborhood  of  t  =  0  at  which  the  first  collapse  occurs, 

d^A  =  -dfzod^A,  (11) 

and  iyAtt  =  c(l  —  ie)“^(V^A  +  |A|*A),  where  c  =  yzQ  >  0.  The  equation 
for  A  now  reads, 

A^  =  +  i\A\’A  -  e(l  -  ic)-^(V2A  +  |A|*A)  .  (12) 

Substituting  the  self-similar  form  (2)  into  (12),  and  using  the  new  vari¬ 
ables  defined  in  (3)  -  (5),  it  is  found  that 

i  5rX  -  X  +  V|x  +  IxI'x  +  /?^^X  = 

1  2 
-([drX  +  +  -y X  +  tx  -  iW  + 

This  equation  gives  the  evolution  of  the  transverse  profile,  Xi  as  long  as  a 
and  /3  are  known. 

For  our  purposes  we  can  obtain  equations  for  a  and  /3  from  the  solvability 
condition  of  the  first-order  solution  in  the  asymptotic  expansion  of  x>  where 

X  =  XO  +  Xi  +  X2  +  • '  •  (  (14) 
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and  it  is  assumed  that  o?  >  drOt,  drX  1)  /?  C  1  as  t  -f  oo.  The 

zero-order  equation  in  this  expansion  is 

v|xo  —  Xo  +  IXol^Xo + P^^Xo  =  0  I  (15) 

and  its  solution  should  satisfy  the  conditions  x{°°)  =  Xf(0)  =  0-  Near  the 
peak  for  small  P(^X0i  (15)  is  approximated  by  the  ground  state  stationary 
solution  of  the  (i-dimensional  NLS  with  7  =  0  in  (1),  Rq. 

The  first-order  equation  is 

LXi  +  ||xo|'"*)(^Xi  = -*«(^)Xo  -  i^rXo 
d  o? 

-  ^[*^(2  -  *(/3  +  — )^^Xo]  >  (16) 

where  ^ 

L=]^v}-i+pe+{^+i)\xov  . 

By  rewriting  (16)  in  terms  of  the  real  and  imaginary  parts  of  xi  using 
=  5  +  iV*,  and  drXo  =  dpxo^rP  one  obtains, 

LjfS  =  -eaip^  -f  id^Xo)  ^  (17) 

=  -9;3Xo5r)3-€[l-(^-|- ^)^^]Xo-u(^)Xo  ,  (18) 


-  1  +  /3e2  +  (1  -H  l)|xo|*  ±  ||xol*-"xE  .  (19) 

It  has  been  assumed  that  f  ^  0(1)  near  p  =  0  and  that  the  term  is 

small.  The  solvability  condition  for  (18)  requires  that 

{Ro\dpXo)drP  =  -v{P)xo  -  e[<i2§)  -  (yS  +  •  (20) 

Using  the  relation,  (iZol^/JXo)  =  (^^•Ro)/‘^»  evolution  of  g  is  given  by  the 
set  of  ordinary  differential  equations 


dr 

^  AO  2 

—  =  4/3  -  , 
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(21) 

(22) 

(23) 


This  reduced  set  of  equations  (21)  -  (23)  describes  the  evolution  of  the 
parameters  of  the  self-similar  solution  near  the  peak  of  the  envelope  in  the 
presence  of  normal  dispersion.  As  in  the  unperturbed  case,  the  behavior  of 
g  indicates  whether  collapse  occurs. 

We  now  analyze  the  phase  portrait  of  (22)  and  (23).  In  Fig.  1  the  phase 
portrait  of  (22)  and  (23)  is  shown  for  c  =  0.1,  s  =  2  and  d  =  2.  In  this  case 
a  =  1.674.  For  any  finite  value  of  the  dispersion  the  point  at  a  =  =  0  is 

no  longer  an  attracting  limit  point.  Trajectories  that  would  have  collapsed 
for  7  =  0  now  tend  to  a  <  0.  Since  g  is  positive  gz  must  be  negative  and 
these  trajectories  must  correspond  to  solutions  that  do  not  collapse,  since 
a  =  gr/g-  At  a  =  0  the  peak  intensity  of  the  pulse  reaches  its  maximum 
and  begins  to  decrease. 

Note  that  this  padr  of  equations  has  two  fixed  points  which  lie  on  the 
parabola  /3  =  o?/A.  The  right  branch  of  this  parabola  has  a  hyperbolic  fixed 
point  while  the  left  has  a  source.  Above  these  points  and  to  the  right  of  the 
left  branch  of  the  parabola  trajectories  tend  to  infinity  along  P  =  o?/A.  It  is 
easy  to  show  that  these  fixed  points  originate  at  infinity  as  c  is  made  finite. 
As  e  is  increased,  this  pair  of  points  moves  toward  the  origin  approaching  a 
point  near  /3  =  2  as  e  tends  to  infinity.  These  fixed  points  are  therefore  not 
relevant  to  the  behavior  of  the  filaments  because  they  are  so  far  away  from 
the  trajectories  along  which  the  initial  surge  towards  collapse  occurs. 

Also  note  that  the  number  of  negatively  dispersive  coordinates,  d,  enters 
into  the  first-order  equation  (16)  through  the  linear  operator  on  the  left  hand 
side  and  only  in  the  second,  almost  purely  real,  term  on  the  right  hand  side. 
Because  of  this,  d,  and  thus  s,  only  enter  the  reduced  equations  through  c 
and  o  and  then  only  in  (23).  Since  the  structure  of  the  phase  portrait  of 
the  reduced  equations  is  relatively  insensitive  to  variations  in  e  and  a  we 
find  the  surprising  result  that  in  addition  to  being  qualitatively  similar,  the 
evolution  of  wave  packets  of  different  dimension  undergoing  critical  collapse 
when  perturbed  by  normal  dispersion  is  also  quantitatively  similar. 

Thus  we  have  shown  analytically  that  a  finite  amount  of  dispersion  has 
the  effect  of  arresting  the  self-similar  collapse  of  a  pulse.  We  have  also 
verified  this  conclusion  using  numerical  simulations  [1]. 

1.2  Threshold  for  collapse 

Although  normal  GVD  can  arrest  critical  collapse,  this  is  not  the  case  if 
the  input  power  is  sufficiently  large,  or  the  GVD  too  small.  To  address  this 
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question  we  return  to  the  NLSE  in  dimensional  form 


(It + Hw) + ^ 


We  have  conducted  extensive  numerical  simulations  of  (24)  with  n2  >  0, 
d^kfOuj"^  >  0,  and  initial  data  of  the  form  A{x,y,z  =  0,t)  =  exp(-(x^  + 
y^)/2wQ  —  t^f2tp),  to  explore  the  dynamics  of  pulse  propagation.  These 
simulations  have  shown  that  as  the  input  power  is  increased  a  threshold, 
Pth>  is  reeiched  where  self  focusing  dominates  the  initial  dynamical  evolution 
of  the  pulse.  Near  PpH  the  pulse  splitting  process  coincides  with  the  arrest 
of  the  collapse  at  the  peak  of  a  pulse. 

The  solutions  of  (24)  are  functions  of  only  two  parameters  when  Gaus¬ 
sian  initial  data  is  used.  We  have  chosen  to  write  the  first  parameter  as 
7  =  Ldf/LdS}  which  manures  the  strength  of  the  dispersion  relative 
to  diffraction,  where  Lds  =  ^l/k"  is  the  linear  dispersion  length,  and 
Ldf  =  Wofc£,/2  is  the  diffraction  length.  The  second  parameter,  p,  mea¬ 
sures  the  strength  of  the  nonlinearity  and  is  the  ratio  of  the  peak  input 
power,  P,  of  a  Gaussian  pulse  to  the  critical  power  for  2D  self  focusing. 
Pc  =  (1.22)^7rA^/(327ion2). 

The  self-focusing  threshold  in  normally  dispersive  media  is  greater  than 
that  for  stationary  2D  self  focusing  so  Pxh  >  Below  PtHi  pulse  splitting 
may  occur  only  after  the  pulse  has  passed  its  focus.  For  P  >  Pth  pulse 
splitting  occurs  at  or  before  the  focus  of  the  pulse.  For  a  series  of  values  of  the 
dispersion  parameter  7,  we  have  numerically  integrated  (24)  to  determine 
the  threshold  power,  Pth-  This  threshold  is  plotted  in  Fig.  2,  where  the 
filled  circles  indicate  simulations  that  were  below  the  threshold,  PxHi  S’lid 
the  open  circles  indicate  simulations  that  were  above. 

An  analytic  estimate  for  the  threshold  power  Pth  be  obtained  us¬ 
ing  physical  arguments.  This  estimate  is  based  on  the  idea  that  when  the 
length  scale  for  self  focusing  is  comparable  to  the  length  scale  for  which 
nonlinear  dispersive  spreading  reduces  the  peak  power  below  the  2D  criti¬ 
cal  power,  self  focusing  should  be  prevented.  The  self-focusing  length  scale, 
zsFf  corresponding  to  the  peak  input  power  in  the  absence  of  GVD  (7  =  0) 

(25) 

This  length  scale  is  an  improved  version  of  that^  conjectured  by  Kelley  on  the 
premise  that  the  collapse  should  occur  near  the  fpcal  point  of  the  nonlinear 
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zsfILdf  ==  0.367 


[(?*- 


0.852)2  _  0.0219 
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lens  induced  by  the  incident  beam.  The  dispersion  length  scale,  znlgvd^ 
taking  into  account  nonlinear  effects,  is-  that  length  over  which  the  combined 
effects  of  self-phase  modulation  and  normal  GVD  would  reduce  the  peak 
power  from  its  input  value  to  the  critictd  power  p  =  1  in  the  absence  of 
diffraction.  To  obtain  zjvlgvDi  we  use  an  approach  previously  developed 
by  Potasek  et.al.  for  pulse  spreading  in  optical  fibers.  A  temporal  analogy 
of  Kelley’s  approach  to  self  focusing  is  employed  whereby  one  considers  the 
“temporal  lens”  induced  by  the  incident  (on-axis)  pulse  which  frequency 
chirps  the  field.  Normal  GVD  then  causes  the  propagating  pulse  width  a{z) 
to  increase  from  its  initial  value  a(0)  according  to 


, 

c{0)  -  \ 


1  +  (^  +  22  )(47^(i:/i/DF)2)  + 


(AjPiz/LDFW-  (26) 


In  passing  from  the  notation  of  Potasek  et.  cJ.  to  ours  we  have  used  (f>r„  = 
2Pz/Ldf^  deff  =  ^yzILDF-  Note  that  when  /3  =  0  the  variation  of  the 
pulse  width  is  due  only  to  linear  spreading.  When  7  is  small  with  respect 
to  /3  the  linear  term  may  be  ignored.  By  assuming  that  the  propagating 
pulse  remains  roughly  Gaussian  so  that  p(0)o’(0)  =  p(z)a(z),  and  setting 
p(0)  =  p,  and  p{z)  =  1,  in  (26),  the  dispersion  length,  znlgvd,  can  be 
written  in  terms  of  p  and  7. 

The  dynamics  of  the  pulse  evolution  are  dictated  by  the  relation  between 
zsF  and  ZffLGVD-  For  ^SF  >  znlgvd  we  expect  normal  GVD  to  lower  the 
peak  power  of  the  pulse  enough  to  prevent  collapse,  sop  <  Pth-  la  contrast, 
for  ZSF  <  zhlgvd  spatial  self  focusing  should  dominate  the  dynamics  and 
will  correspond  to  p  >  Pth-  -A-n  estimate  for  the  threshold  power  Pth 
in  terms  of  the  dispersion  parameter  7  can  be  obtained  by  letting  zsf  = 
znlgvd-  Excluding  the  linear  contributions. 


znlgvd/Ldf  = 


[3.38-h5.2(p2-  1)]*  -1.84 


157P 


and. 


(27) 


I  [3.38  -I-  5.2{P^jj  -  1)]  *  -  1.84|  [(P/h  -  0-852)^  -  0.0219] 


2Pth 


(28) 


This  analirtic  estimate  is  displayed  as  the  solid  line  in  Fig.2.  It  provides  an 
excellent  characterization  of  the  numerical  data  ipver  a  wide  range  of  7.  It 
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works  well  for  both  7  <  1,  which  is  the  experimentally  important  parameter 
regime,  and  for  7  ~  1.  If  the  linear  spreading  distance,  zqvDi  is  used,  where 
=  0  in  (26),  one  obtains  the  curve  annotated  with  dots  and  dcishes,  ( — ), 
plotted  in  Fig.  2. 

1.3  Other  issues 

The  main  findings  of  this  section  are  that  normal  GVD  can  arrest  collapse 
for  powers  close  to  critical,  but  only  for  powers  below  a  threshold  value 
dependent  upon  the  value  of  the  GVD.  Below  this  threshold  the  collapse  is 
arrested  via  pulse-splitting,  but  above  the  threshold  even  the  pulse  splitting 
is  not  sufficient  to  arrest  collapse  before  the  envelope  approximations  fail. 

In  view  of  the  importance  of  these  results  we  have  investgated  the  col- 
laspe  further  to  look  for  potential  experimental  signals  of  the  pulse-splitting 
phenomenon  [2].  Our  numerical  simulations  show  that  conical  emission  ac- 
companys  pulse-splitting,  the  frequency  of  the  emitted  cone  being  depen¬ 
dent  on  the  value  of  the  GVD.  The  conical  emission  is  also  accompanied  by 
significant  superbroadening  of  the  pulse,  and  is  responsible  for  the  pulse¬ 
splitting  phenomenon.  We  have  related  both  the  spectral  super  broadening 
and  conical  emission  to  a  fundamental  four-wave  interaction  which  can  be¬ 
come  phase-matched  via  the  normal  GVD. 

We  have  also  begun  to  explore  the  arrest  of  SF  collapse  beyond  the  en¬ 
velope  approximation.  For  this  we  have  developed  Maxwell  equation  solver 
which  are  capable  of  handling  nonlinearity.  Numerical  simulations  indicate 
that  vectorial  field  effects  arrest  collapse  as  the  field  approaches  wavelength 
scales,  even  in  the  absence  of  normal  GVD.  However,  there  is  no  support¬ 
ing  theory  for  this  at  present.  We  have  also  discovered  a  new  phenomenon 
which  is  beyond  the  envelope  approximation,  namely,  carrier  wave  shocking 
[4].  Optical  shocks  have  been  studied  for  sometime  in  the  area  of  nonlinear 
fiber  optics.  These  shocks  appear,  however,  in  the  field  envelope.  We  have 
found  that  for  ultrahort  pulses  shocks  can  start  to  appear  on  the  carrier, 
and  is  characterized  by  the  appearence  of  a  strong  third-harmonic.  Carrier 
shocks  may  play  a  significant  role  in  the  arrest  of  SF  collapse  when  the  pulse 
length  contracts  and  the  intensity  gets  very  large,  since  the  energy  converted 
to  the  thrid-harmonic  may  be  sufficient  to  reduce  the  fundamental  power 
below  critical. 
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2  Self-focusing  and  Optical  Breakdown 

In  this  section  we  describe  our  achievments  in  the  area  of  modelling  of  SF 
and  LIB  in  water  as  a  model  for  the  vitreous  humor.  The  full  details  of  our 
extensive  numerical  studies  are  presented  in  Refs.  [5,  6],  and  here  we  display 
the  most  revealing  results.  Research  highlights  include 

•  Development  of  a  comprehensive  model  for  nonlinear  propagtion  in¬ 
cluding  SF,  GVD,  LIB,  and  MPA  [6]. 

•  Development  of  a  code  for  nonlinear  propagation  which  was  transferred 
to  Brooks  Air  Force  Base  for  studies  of  eye  damage. 

•  Analytic  criteria  for  the  relative  roles  of  SF  and  LIB  as  a  function  of 
pulse  length  and  focusing  conditions  [5]. 

•  Calculations  of  breakdown  thresholds  with  and  without  SF  and  com¬ 
parison  with  previous  theories  and  experiments  where  possible.  This 
provided  validation  of  the  model. 

•  Demonstration  that  SF  and  LIB  are  both  key  players  in  ultrashort 
pulse  propagtion  in  water. 

•  Detailed  investigation  of  strongly  nonlinear  behavior  of  ultrshort  pulse 
collapse  shwoing  pulse  splitting  and  potentially  new  mechanisms  for 
eye  damage. 

2.1  TheoreticeJ  Model 

To  start  we  describe  the  general  theoretical  model  for  for  nonlinear  propa¬ 
gation  in  water  which  is  a  generalization  of  the  model  previously  studied  by 
Feit  and  Fleck  [10]  to  include  GVD,  SF,  and  MPA.  Then  assuming  propa¬ 
gation  along  the  z-axis,  the  equation  for  the  electric  field  envelope  5(r,  z,  t) 
in  a  reference  frame  moving  at  the  group  velocity  is 

where  the  terms  on  the  right-hand-side  describe  transverse  beam  diffraction, 
GVD,  nonlinear  SF,  absorption  and  defocusin^  due  to  the  electron  plasma, 
and  MPA  involving  K  photons,  respectively.  Here  w  is  the  optical  frequency, 
\S\'^  the  intensity,  k  =  n^fco  =  nbu/c,  the  quanti!j;y  k"  =  d’^kfdu^  controls 
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the  magnitude  and  sign  of  the  GVD,  with  k"  >  0  corresponding  to  normal 
dispersion  and  k"  <  0  to  anomolous  dispersion,  712  is  the  nonlinear  coef¬ 
ficient  such  that  the  nonlinear  change  in  refractive-index  is  n2|5p,  a  the 
cross-section  for  electron-neutral  inverse  bremsstrahlung,  r  is  the  electron 
collision  relaxation  time,  and  is  the  nonlinear  coefficient  for  if-photon 
absorption.  The  order  of  the  MPA  is  obtained  from  K  =  mod{Eg/fujj)  +  1, 
which  is  the  minimum  number  of  photons  of  energy  tiuj  needed  to  overcome 
the  ionization  energy  Eg  for  liberating  an  electron  from  a  water  molecule. 

The  Drude  model  equation  which  describes  the  electron  density  p  is  [11] 


dt 


1  a 


nj  Eg 


Khu) 


—  ap 


(30) 


The  first  term  on  the  right-hand-side  of  this  equation  describes  growth  of  the 
electron  plasma  by  cascade  (avalanche)  ionization,  the  second  term  is  the 
contribution  of  MPA,  which  acts  both  as  a  source  for  the  cascade  process  and 
as  a  contributor  to  plasma  growth,  and  the  third  term  describes  the  radiative 
electron  recombination.  In  addition,  if  there  is  a  free-electron  density  in  the 
water  prior  to  the  pulse  it  will  also  be  amplified  by  cascade  ionization.  Thus 
we  also  include  a  background  initial  density  ph  as  an  initial  condition  for  Eq. 
(30). 

Fig.  3(a)  shows  the  focusing  geometry  we  considered.  We  are  interested 
in  the  solutions  of  Eqs.  (29)  and  (30)  for  an  initial  collimated  Gaussian 
beam  which  enters  the  water  sample  following  a  lens  of  focal  length  / 


5(.,0,t)=gexp(-^-|-i^),  (31) 


where  Pin  is  the  peak  input  power  [P(i)  =  j2Trrdr\£(r,z  =  0,  = 

P»nexp(— 2t^/tp)],  Wo  the  spot  size,  and  Tp  =  2fp  is  the  full  temporal  width 
at  the  1/e^  points  of  the  pulse  intensity  distribution. 

We  have  considered  a  specific  focusing  geometry  with  /  =  1.7  cm,  this 
being  roughly  the  focal  length  of  the  eye.  Unless  otherwise  stated  the  input 
spot  size  was  taken  as  wq  =  200  pm,  which  focuses  at  d  =  1.69  cm  with 
Wf  =  11.8  pm.  The  Rayleigh  range  of  the  focused  beam  is  zy  =  0.1  cm.  The 
remaining  parameters  used  for  water  for  a  wavelength  of  Aq  =  580  nm  are 
Ub  =  1.33,  712  =  4.1  X  10"^®  cmVW,  K  =  A,  =  1.55  x  lO"**®  m^W"®, 
Eg  =  6.5  eV,  T  =  10~^®  s,  iiC  =  4,  a  =  1.4  X  10“^^  cm®,  k"  =  0.05  ps®/m, 
and  Per  =  IMW. 
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2.2  Light-induced  breakdown  with  linear  focusing 

To  provide  a  framework  in  which  to  view  the  numerical  simulations  with 
SF  and  LIB,  we  first  looked  at  LIB  with  no  SF.  This  is  closely  realized 
experimentally  by  the  focusing  geometry  shown  in  Fig.  3(b).  In  the  absence 
of  SF  we  may  simply  solve  Eq.  (30)  for  the  density  at  the  linear  focus  of 
the  beam 

p{t)  =■  pbgip)  d"  Pmpip)}  (32) 

where 


Pbg{t)  =  pb<p(t),  Pmp{t)  =  <P(t)^fL  du, 


<p(t)  =  exp 


IL 


Pbg{t)  is  the  homogeneous  part  of  the  solution  which  arises  from  the  back¬ 
ground  electron  density  pi,  and  pmpify  is  the  inhomogeneous  solution  which 
grows  from  the  electron  density  initiated  by  MPA.  Since  the  analytic  so¬ 
lution  (32)  is  not  very  informative,  we  have  numerically  solved  the  density 
equation  (30)  to  obtain  the  threshold  for  LIB.  Here  we  use  as  our  criterion 
for  LIB  that  p{t  =  oo)  =  pbd  =  10^®  cm~®,  where  p{t  =  oo)  is  the  final 
density  following  the  pulse.  This  criterion  for  LIB  was  selected  to  coincide 
with  the  appearence  of  vapour  bubbles  in  the  water  sample,  also  known 
as  the  bubble  endpoint.  The  breakdown  intensity  versus  pulse  duration  is 
shown  in  Fig.  4(a),  and  the  breakdown  power  in  Fig.  4(b),  for  the  three 
backround  densities  pt  =  0  cm“^  (dotted  line),  10^®  cm~^  (solid  line),  and 
10^^  cm“®  (dashed  line).  The  threshold  intensity  and  power  are  seen  to  be 
insensitive  to  the  background  density  for  pulse  durations  less  than  10  ps, 
and  the  variation  in  the  predicted  threshold  is  only  a  factor  of  two  even  for 
pulses  as  long  as  200  ps. 

Table  I  shows  a  comparison  of  the  breakdown  intensity  Ibd  calculated 
here  (fourth  column)  with  the  previous  experimental  (second  column)  and 
first-order  theoretical  (third  column)  results  of  Kennedy  et.  al.  [12,  13], 
and  Hammer  et.  al.  [14],  with  good  overall  agreement.  The  experiment  of 
Hammer  et.  al.  [14]  was  performed  using  a  focusing  geometry  as  in  Fig. 
3(b)  with  a  thin  water  sample  and  at  a  wavelength  of  580  nm.  Comparison 
with  the  experimental  data  provides  verification  of  the  first-order  model 
and  also  of  the  results  from  the  Drude  model  contained  in  EJq.  (30).  This 
comparison  gives  validation  of  the  Drude  model  as  a  reasonable  model  for 
ultrashort  pulse  LIB.  v 
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2.3  Multiphoton  initiation  versus  background  electron  den¬ 
sity 

The  calculations  perfornaed  include  components  of  the  generated  electron 
density  which  are  initiated  both  from  the  background  electron  density  pb 
and  from  MPA.  Both  of  these  seed  electron  sources  are  amplified  by  cascade 
ionization.  To  determine  which  of  these  is  the  dominant  seed  electron  source 
we  have  calculated  the  ratio  of  the  background  contribution  to  the  total 

C 

By  evaluating  this  ratio  as  a  function  of  pulse  duration  we  find  that  C  =  1/2 
for  Tp  =  Tx  =  160  ps,  so  that  MPA  and  the  backround  components  are  in 
exact  balance.  This  value  is  further  supported  by  the  first-order  model  [12] 
which  3rields  a  pulse  duration  of  Tx  =  120  ps  for  the  cross-over.  For  pulse 
durations  less  than  the  cross-over,  Tp  <  Tx,  MPA  dominates  the  cascade 
initiation,  and  for  Tp  >  Tx,  the  background  density  dominates  the  cascade 
initiation. 

2.4  Numerical  simulations  with  self-focusing 

The  simulations  were  performed  using  the  beam  propagation  method  in 
which  the  propagation  is  broken  into  small  linear  and  nonlinear  steps  which 
are  performed  sequentially.  Considerable  advantage  is  gained  by  making 
use  of  the  fact  that  for  P  <  Per  or  P  ~  Per.  in  the  initial  stage  of  the 
propagation,  when  the  nonlinear  effect  is  not  important  and  the  plasma  has 
not  been  generated,  the  propagation  is  described  by  the  linear  equation  with 
only  diffraction  and  dispersion.  The  solution  of  this  equation  can  be  written 
down  exactly, 

5(r,  f,  z)  =  SoDfir,  z)D.(t,  z),  (34) 

where 


We  use  this  linear  equation  to  propagate  the  initial  field  into  the  nonlinear 
interaction  zone,  which  is  then  used  as  the  input  tp  the  nonlinear  code.  This 


Pbgjoo) 

Pbgioo)  +  Pmpioo)' 
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represents  a  considerable  reduction  in  computation  time,  and  also  allows 
us  to  use  the  computational  grid  more  effectively,  since  the  bulk  of  the 
linear  focusing  by  the  lens  is  dealt  with  before  passing  over  to  the  nonlinear 
propagation  code. 

Turning  now  to  the  numerical  simulations  with  SF,  we  have  determined 
the  threshold  for  breakdown  in  the  same  manner  as  before.  The  results  of 
our  numerical  simulations  are  summarized  in  Fig.  5  where  we  show  (a)  the 
breakdown  threshold  of  peak  input  power  {Pin),  normalized  to  the  critical 
power  for  collapse,  and  (b)  the  threshold  intensity,  both  as  functions  of  pulse 
duration  Tp.  The  solid  circles  in  Fig.  5(a)  are  our  data  points  corresponding 
to  a  peak  plasma  density  of  10^®  cm“®.  The  solid  curve  is  the  threshold 
power  curve  obtained  in  the  previous  Section  without  SF,  and  also  shown 
in  Fig.  4(b).  We  see,  therefore,  that  for  pulses  longer  than  Tp  =  10  ps, 
the  thresholds  with  and  without  SF  agree  very  well.  This  is  to  be  expected 
since  for  these  pulse  durations  the  breakdown  power  is  considerably  below 
the  critical  power,  Pbd  <  Per,  and  SF  plays  only  a  minor  role.  In  contrast, 
for  Tp  <  10  ps  the  breakdown  power  is  seen  to  tend  towards  the  critical 
power  Per.  This  is  to  be  expected  intuitively,  since  as  the  input  power 
approaches  the  critical  power  small  changes  in  the  input  power  can  cause 
massive  changes  in  the  peak  intensity  at  the  laser  beam  focus,  which  in  turn 
can  cause  very  large  changes  in  the  peak  plasma  density  generated,  so  that 
the  plasma  density  for  breakdown  is  easily  reached.  This  is  particularly  true 
of  the  component  of  the  density  initiated  by  MPA  since  the  source  term 
varies  as  \S\^^  (see  Eq.  (30)),  with  K  =  4,  and  MPA  is  expected  to  be  the 
dominant  initiation  mechanism  for  plasma  generation  for  pulse  durations 
less  than  Tx  =  160  ps. 

The  corresponding  breakdown  intensities,  chosen  as  the  peak  intensities 
which  appeared  during  the  respective  numerical  simulations,  versus  pulse 
duration  are  shown  in  Fig.  5(b)  (solid  circles)  along  with  the  result  with¬ 
out  SF  from  Fig.  4(a)  (solid  line).  Here  we  see  a  remarkable  agreement, 
especially  for  longer  pulses,  between  the  results  with  and  without  SF.  This 
agreement  stems  from  the  fact  that  LIB  is  triggered  by  the  local  electric 
field  strength,  or  intensity,  so  that  LIB  has  an  intensity  threshold  [11]  which 
is  essentially  independent  of  the  focusing  conditions. 

Plasma  shielding  arises  when  the  incident  pulse  energy  is  absorbed  and 
converted  to  the  generation  of  the  electron  plasma  via  LIB,  thus  shielding 
the  region  beyond  the  plasma  from  the  full  pulse  energy.  In  Fig.  6  we  show 
the  peak  generated  plasma  density  and  corresponding  percentage  energy 
absorbed  for  pulse  durations  Tp  =  200  ps  (solid  line),  100  ps  (dotted  line). 
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and  20  ps  (dashed  line).  In  each  case  the  lowest  power  indicated  is  slightly 
above  the  threshold  power  for  LIB,  with  weak  plasma  shielding  in  each 
case.  However,  far  above  the  threshold  the  peak  plasma  density  and  the 
plasma  shielding  both  increase,  with  10-25  %  absorption  depending  on  the 
pulse  length.  This  clearly  demonstrates  that  plasma  shielding  becomes  more 
significant  above  the  bubble  endpoint  threshold,  as  expected  on  physical 
grounds. 

An  interesting  aspect  of  our  model  is  the  behavior  of  nonlinear  prop¬ 
agation.  Before  the  high-density  plasma  is  generated,  the  propagation  is 
essentially  described  by  the  self-focusing  solution  However,  once  the  plasma 
it  acts  to  absorb  energy  and  defocus  the  beam.  The  process  of  absorption 
and  defocusing  could  result  in  very  rich  nonlinear  phenomena  which  have 
not  been  studied  in  detail.  Here  we  report  a  case  of  numerical  simulation 
that  shows  qucditatively  different  behavior  from  what  we  have  described  in 
the  previous  section.  We  chose  a  smaller  input  spot  size,  wq  =  120  /xm,  so 
as  to  reduce  the  effect  of  linear  focusing.  This  allows  us  to  follow  the  nonlin¬ 
ear  self-focusing,  which  is  very  explosive,  within  the  limit  of  our  code.  The 
process  of  pulse  propagation  and  plasma  generation  is  presented  in  Figs.  7 
and  8.  Fig.  7  shows  maa;t|e(r  =  0,  z,  t)/fio|  (top),  maxtp{r  =  0,  z,  t)  (middle) 
and  energy  (bottom)  as  functions  of  propagation  distance  z.  It  is  interesting 
to  note  that  there  are  two  peaks  in  maxtp{r  =  0,z,t),  one  in  front  of  and 
the  other  behind  the  focal  point  z  =  f.  The  spatial  extension  of  the  plasma 
generated  is  wider  than  all  the  cases  described  in  previous  sections.  One  can 
clearly  see  that  the  energy  absorption  occurs  in  a  region  where  the  plasma  is 
generated,  while  in  other  regions  of  low  plasma  density,  the  energy  is  almost 
constant.  Fig.  8  shows  p  (top)  and  e  (bottom)  as  functions  of  (x,  t)  at  three 
propagation  distances:  (a)  z  =  1.6869,  (b)  z  =  1.70043,  and  (c)  z  =  1.73562. 
The  pulse  starts  to  split  in  front  of  the  focal  point.  While  the  trailing  split 
piece  (the  shorter  piece  corresponding  to  higher  plasma  density)  is  damped 
once  the  splitting  occurs,  the  leading  piece  continues  to  focus  and  generates 
high  plasma  density  behind  the  focal  point  [see  Fig.  8(c)].  Although  here 
the  energy  absorption  is  relatively  higher,  the  propagating  pulse  shows  a 
tendency  of  not  being  absorbed.  This  suggests  a  possible  damaging  mecha¬ 
nism:  The  leading  split  piece  of  the  pulse  may  escape  the  plasma  generated 
in  front  of  the  focal  point  and  reach  the  retina,  causing  additional  damage. 
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3  Nonlinear  Absorbing  Interface 

This  section  describes  our  work  on  nonlinear  absorbing  interfaces  as  a  start 
to  addressing  the  issue  of  nonlinear  effects  at  the  boundary  between  the 
transparent  vitreous  humor  and  the  highly  absorbing  retinal  pigment  ep¬ 
ithelium.  The  models  investigated  are  not  directly  relevant  to  the  eye  but 
shine  a  great  deal  of  light  on  the  issues  which  arise  at  highly  absorbing  inter¬ 
faces  when  the  absorption  per  wavelength  approaches  unity.  The  key  point 
is  that  in  this  regime  the  envelope  approximation  is  invalid  and  an  approach 
based  on  Maxwell’s  equations  must  be  employed.  Previous  work  in  this  area 
concentrated  only  on  the  steady-state  properties  whereas  our  interest  was 
in  ultrshort  pulses.  Research  highlights  include 

•  Theoretical  prediction  of  a  new  effect  whereby  an  incident  pulse  is 
Doppler-shifted  as  a  result  of  reflection  of  a  moving  absorption  front 
[8] 

•  Demonstration  of  induced  focsuing  of  an  ultrashort  pulse  reflected  of 
a  saturable  absorber  [9]. 

Here  we  concentrate  solely  on  the  new  results  on  the  nonlinear  skin  effect. 
3.1  Nonlinear  skin  effect 

In  the  linear  optical  skin  effect  a  pulse  incident  from  air  is  reflected  from  a 
highly  absorbing  interface  after  penetrating  only  a  fraction  of  a  wavelength 
into  the  absorbing  medium,  this  distance  being  the  skin  depth  The  skin 
effect  is  therefore  of  fundamental  importance  in  understanding  the  electro¬ 
dynamics  of  pulse  propagation  at  condensed  matter  interfaces.  In  addition, 
it  belongs  to  an  important  class  of  optical  problems  for  which  the  notion  of 
an  electromagnetic  field  envelope  varying  slowly  on  the  scale  of  a  wavelength 
simply  does  not  apply. 

We  have  discovered  a  nonlinear  generalization  of  the  skin  effect.  Broadly 
speaking,  saturation  of  the  absorption  allows  the  incident  field  to  penetrate 
beyond  the  linear  skin  depth  into  the  medium,  and  this  causes  an  absorp>- 
tion  front  to  propagate  into  the  medium  which  separates  the  regions  of  low 
(saturated)  and  high  (unsaturated)  absorption.  The  front  is  excited  by  the 
incident  pulse  which  is  in  turn  reflected  fron^  the  sharp  absorption  front, 
yielding  a  self-reflected  pulse.  Thus  the  absorption  front  acts  as  a  moving 
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mirror  from  which  the  pulse  is  self-reflected,  and  the  pulse  suifers  a  red-shift 
due  to  the  Doppler  effect. 

We  consider  the  time-dependent  propagation  of  a  linearly  polarized  plane 
electromagnetic  wave  incident  on  a  nonlinear  medium  composed  of  two-level 
systems.  For  propagation  along  the  z-axis,  and  taking  the  electric  field 
polarized  along  the  x-axis,  Maxwell’s  curl  equations  take  the  form 

dBy  dE^  dD^  dHy  . 

dt  ~  dz'  dt  ~  dz'  ^  ^ 

where  By  =  fioHy.  The  specification  of  the  problem  is  completed  with  the 
constitutive  relation  =  cqEx  +  Px,  where  Px  is  the  optical  polarisation. 
To  elucidate  the  basic  physics  we  employ  a  two-level  model  to  describe  the 
optical  response  with  lower  electronic  state  |1)  and  upper  state  |2).  The 
Bloch  equations  are  then 

-^  +  (T2  +  «*'2i)/>21  =  +  =  2i^(/>2i  —  P2i)>  (36) 

where,  p2\  is  the  off-diagonal  density  matrix  element,  n  =  (pn  —  P22)  is 
the  population  difference  between  the  lower  and  upper  states,  W21  is  the 
transition  frequency,  p  is  the  dipole  moment  in  the  field  direction,  and  71 
and  72  are  phenomenological  damping  constants  for  the  population  and 
polarisation,  respectively.  The  polarization  due  to  the  atoms  is  then  given 
by  Px  =  N{z)p{p2i-\-c.c.),  with  N{z)  the  density  of  two-level  systems  which 
varies  along  z  in  general. 

We  have  solved  the  coupled  Maxwell-Bloch  Eqs.  (35)  and  (36)  using  a 
standard  discretization  of  the  Maxwell  equations  as  described  by.  The  Bloch 
equations  were  integrated  in  time  using  a  standard,  fourth-order  Runge- 
Kutta  method.  The  system  of  equations  was  solved  with  the  initial  condition 
on  the  field 

Ex{z,  t  =  0)  =  £7o  cos(27ra;(z  -  Z())/c)e-(*-'“)*/("‘'>>*,  (37) 

along  with  a  similar  expression  for  By  with  Bq  =  Eq/ c.  Here  Eq  is  the  peak 
input  electric  field,  w  =  27rc/Ao  is  the  central  pulse  frequency,  tq  =  2to  is 
the  full-width  at  the  1/c^  points  of  the  pulse  intensity  profile  in  time  units, 
and  zq  is  the  position  of  the  pulse  center  at  t  =  0.  The  nonlinear  interface 
Wcis  imposed  by  tailoring  the  density  profile  N{z)  =  d{z  —  Ziru)No,  with  z^nt 
the  longitudinal  position  of  the  interface.  The, medium  was  initialized  using 
P2i  =  0)  S'lid  n  =  1  in  the  medium.  When  initializing  the  field  we  ensured 
that  the  field  protruded  negligibly  into  the  nonlinear  medium  at  t  =  0. 
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In  the  limit  of  cw  fields,  as  previously  studied  by  Roso-Franco,  the  self- 
reflected  wave  arises  when  the  normalized  parameters, 


-  p_  pEp 


(38) 


are  both  greater  than  unity.  Physically,  V*  determines  the  linear  absorption 
per  wavelength,  oAq  =  27rV’,  and  this  quantity  should  be  greater  than  unity 
for  the  linear  skin  effect.  For  saturation  of  the  absorption  and  self-reflection 
F  >  1  is  required,  since  F^  is  the  peak  incident  field  intensity  normalized 
to  the  cw  saturation  intensity. 

We  consider  the  transient  regime  in  which  the  incident  pulse  width  to  is 
much  shorter  than  the  population  relaxation  time  Ti  =  7f  but  longer  than 
the  polarization  dephasing  time  T2  =  7^^,  Ti  »  to  >  T2.  For  concreteness 
we  adopt  the  following  specific  parameters,  w  =  W12  =  2  X  10^®  rads“^, 
to  =  300  fs,  Ti  =  0.5  ns,  T2  =  50  fs,  p  =  5eao  =  4  x  10-^®  Cm,  JNTo  =  4  X 10^® 
cm“®  and  zq  =  —225  pm.  For  these  parameters  V’  =  3.8  so  that  the  linear 
skin  effect  is  expected  at  low  input  intensities.  We  have  numerically  verified 
that  this  is  indeed  the  case,  and  the  input  pulse  suffers  minimal  distortion 
in  profile  or  spectrum  upon  reflection. 

Figures  9(a)  and  1(b)  show  an  example  of  the  calculated  pulses  at  two 
different  times  for  a  peak  input  field  of  Eq  =  1.6  X  10®  V/m.  Although 
our  calculations  employ  the  full  field,  we  display  only  the  envelope  obtained 
from  joining  the  peaks  as  shown  by  the  solid  line,  since  it  is  not  possible 
to  resolve  the  carrier  in  the  plots.  The  field  strength  is  associated  with  the 
sccde  shown  on  the  left-hand-side  of  the  plots.  In  Fig.  9(a)  for  t=1.68  ps 
the  peak  of  the  input  pulse  has  not  yet  reached  the  interface  at  Zint  =  400 
pm,  but  one  can  clearly  see  the  leading  edge  of  the  pulse  is  penetrating 
only  a  short  distance  into  the  interface,  as  expected  for  the  skin  eflfect.  The 
dashed  line  in  Fig.  9(a),  which  is  associated  with  the  right-hand  scale,  is  the 
local  wavelength  for  the  field.  This  is  determined  numerically  by  calculating 
the  local  wavenumber  K  via  w  —E'^/Ex,  where  a  prime  signifies  a  z- 
derivative,  and  converting  to  wavelength.  In  Fig.  9(a)  the  local  wavelength 
remains  constant  at  the  input  value  Aq  =  942  nm.  In  contrast.  Fig.  9(b) 
shows  the  field  profile  at  a  later  time  t  =  2.75  ps  following  reflection  from 
the  interface  (for  times  between  the  results  shown  in  Figs.  9(a)  and  9(b),  the 
field  profile  shows  strong  ringing  due  to  interference  between  the  incident  and 
reflected  fields).  The  reflected  pulse  has  developed  a  double-peaked  structure 
(solid  line),  and  become  significantly  chirped ''(dashed  line).  In  particular, 
the  central  portion  of  the  pulse  has  a  peak  local  wavelength  of  990  nm. 
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a  significant  red-shift.  This  red-shift  is  also  evident  in  the  refiected  pulse 
spectrum  (solid  line)  shown  in  Fig.  9(c)  corresponding  to  Fig.  9(b),  along 
with  considerable  spectral  broadening  and  modulation  (the  input  spectrum 
is  shown  by  the  dashed  line  and  is  eissociated  with  the  left-hand  scale). 

The  results  shown  in  Fig.  9  are  typical  of  what  we  observe  in  our  sim¬ 
ulations  in  the  nonlinear  regime,  namely,  distortion  of  the  reflected  field 
profile  and  significant  spectral  modulation  and  associated  red-shift.  To  ex¬ 
pose  the  physics  underlying  these  phenomena  we  show  in  Fig.  10(a)  and 
10(b)  the  spatial  distribution  of  the  full  field  and  the  population  difference 
n  =  —  Pn)  at  various  times  corresponding  to  the  results  shown  in  Fig. 

9.  Figure  10(a)  shows  that  the  field  penetrates  progressively  further  into 
the  interface  as  the  absorption  is  saturated,  as  can  be  seen  by  comparing 
the  field  profiles  at  t  =  1.61  ps  (solid  line)  and  t  =  1.95  ps  (dashed  line)  or 
t  =  2.28  ps  (single-dash-dot  line)  (at  t  =  2.62  ps  the  field  has  been  mostly 
refelcted).  At  t  =  1.61  ps  (solid  line)  Fig.  10(b)  shows  that  n  =  0  before 
the  interface  signifying  zero  absorption,  and  n  =  1  beyond  the  interface  sig¬ 
nifying  large  absorption  due  to  the  two-level  systems.  At  later  times  after 
the  input  pulse  has  penetrated  into  the  interface,  the  population  difference 
is  depleted  and  the  absorption  front  is  seen  to  propagate  into  the  nonlinear 
medium.  The  propagating  absorption  front  maintains  a  sharp  wavelength 
scale  transition  region  so  that  the  linear  skin  effect  still  occurs  but  now 
from  a  moving  absorption  front.  Thus  the  self-reflected  field  must  suffer  a 
red-shift  due  to  the  Doppler-effect,  akin  to  reflection  from  a  mirror  moving 
away  from  a  source.  To  vdidate  this  physical  picture  we  have  determined 
the  absorption  front  velocity  from  the  numerical  simulation  in  Fig.  10(b), 
and  the  result  is  shown  in  Fig.  10(c).  After  initially  accelerating  the  front 
reaches  a  maximum  velocity  of  Vmoi/c  =  0.023  before  decelerating  back  to 
zero  velocity.  The  maximum  wavelength  shift  of  the  reflected  pulse  due  to 
the  Doppler  effect  is  then  AA/Aq  =  2v,„aa,/c  for  Umo*/ c  -C  1  or  AA  =  43  nm 
for  the  free-space  wavelength  Aq  =  942  nm  used  here.  Thus,  based  on  the 
Doppler  eflfect  we  expect  a  maximum  local  wavelength  of  A  =  985  nm,  in 
good  agreement  with  the  numerical  results  in  Fig.  10(b).  The  Doppler  effect 
upon  reflection  from  the  moving  absorption  front  can  therefore  explain  the 
magnitude  of  the  observed  pulse  wavelength  chirp. 

To  summarize,  we  have  introduced  the  dynamic  nonlinear  optical  skin 
effect  for  reflection  of  pulses  from  a  highly  absorbing  interface.  This  new 
basic  effect  for  the  electrodynamics  of  interfaces  combines  the  concepts  of 
self-reflected  waves  and  front  propagation,  anti  is  also  a  prime  example  of  a 
nonlinear  optical  phenomenon  where  the  SVEAffails  and  the  full  Maxwell 
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equations  must  be  employed.  We  have  shown  that  the  nonlinear  optical 
skin  effect  arises  from  moving  absorption  fronts  so  that  the  red-shifting  and 
spectral  modulation  of  the  reflected  pulse  are  clear  experimental  signatures 
of  the  effect. 
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5  Interactions/Presentations 

5.1  Pzirticipation/Presentations 

•  “Nonlinear  Mechanisms  for  Laser  Induced  Breakdown  of  Femtosecond 
Pulses",  3rd  AFOSR  Ultrashort  Laser  Program  Collaborative  Work¬ 
shop,  Brooks  Air  Force  Base,  San  Antonio,  Texas.  14-16  December 
1994. 

•  “Issues  Relating  to  the  Arrest  of  Critical  Collapse  of  Femtosecond 
Pulses  in  Normally  Dispersive  Self-Focusing  Media”,  invited  talk  - 
Conference  on  Nonlinear  Coherent  Structures  in  Physics  and  Biology, 
10-14  July  1995,  Heriot-Watt  University.  Edinburgh. 
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•  “Arrest  of  Critical  Collapse  of  Ultra-Short  Pulses  in  Optically  Trans¬ 
parent  Media”,  Physics  Colloquium,  Georgia  Tech.,  24  April  1995. 

•  “Issues  in  the  Arrest  of  Critical  Collapse  of  Femtosecond  Optical  Pulses 
in  Optically  Transparent  Media”,  Physics  Colloquium,  Free  University, 
Amsterdam,  28  June  1995. 

•  “Computational  Nonlinear  Optics” ,  Computational  and  Physical  Math¬ 
ematics  Contractors  Meeting,  Phillips  Laboartory,  Kirtland  AFB,  28- 
30  June  1995. 

5.2  Consultative  and  Advisory  Functions  to  Other  Labora¬ 
tories  and  Agencies 

•  Visit  by  J.  Moloney  (UA)  to  “3rd  AFOSR  Ultrashort  Laser  Program 
Collaboartive  Workshop”,  Brooks  Air  Force  Base,  San  Antonio,  Texas, 
14-16  December  1994,  to  consult  with  Dr  R.  Albanese  (AL/OES), 
Dr  M.  Potasek  (AL/OES),  Dr  M.  Rogers  (AL/CA),  Dr  B.  Rockwell 
(AL/OEOP),  Dr  P.  Kennedy  (AL/OEOP),  Dr  C.  Cain  (AL/OEOP) 
and  Dr  R  Thomson  (OpTech). 

•  Brooks  Air  Force  visit  to  ACMS,  Department  of  Mathematics,  Uni¬ 
versity  of  Arizona,  9-12  May  1995.  Personnel  involved  included:  J. 
Moloney  (UA),  E.  Wright  (UA),  Q.  Feng  (UA),  R.  Thompson  (Op 
Tech),  P.  Kennedy  (AL/OEO)  and  M.  Potasek  (AL/OES). 

•  Brooks  Air  Force  visit  to  ACMS,  Department  of  Mathematics,  Univer¬ 
sity  of  Arizona.  17-20  August  1995.  Discussion  on  modeling  nonlinear 
propagation  of  picosecond  and  femtosecond  laser  pulses  in  water  and 
in  the  eye.  Personnel  involved  included:  J.  Moloney  (UA),  E.  Wright 
(UA),  Q.  Feng  (UA),  R.  Thompson  (Op  Tech),  P.  Kennedy  (AL/OEO) 
and  M.  Potasek  (AL/OES). 

•  Hosting  and  organisation  of  the  AFOSR/ACMS  “Nonlinear  Optics” 
Workshop,  Tucson,  Arizona.  9-11  October  1995. 


This  workshop  acts  a  a  contract  review  for  AFOSR  and  brings  Air 
Force  Laboartory  scienctists  in  direct  contact  with  academic  researchers 

•  Deliverable:  Computer  Code  to  simulate  critical  self-focusing  of  pi¬ 
cosecond  pulses  in  ocular  and  related  med(a  with  plasma  generation 
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through  avalanche  ionization  included.  Code  is  currently  being  run  by 
Dr  Randy  Thompson  (Op  Tech)  and  Dr  Paul  Kennedy  (AL/OEO)  at 
the  Armstrong  Laboratory. 
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Table  I:  Breakdown  threshold  intensities  (W/cm^). 


TP 

Ibd^ 

2.4  ps 
400  fs 
100  fs 

5  X  10“ 
1.26  X  10^2 
5.6  X  10^2 

6.66  X  10“ 
2.56  X  10^2 
5.46  X  10^2 

1.  Experimental  thresholds  for  50%  probability  of  breakdown 

2.  First-order  calculation  of  the  intensity  breakdown  threshold. 

3.  Intensity  breakdown  threshold  from  present  work. 


FIGURE  CAPTIONS 

Fig.  1  The  phase  portrait  of  the  reduced  equations  (22)  and  (23)  for  s  =  2 
and  d  =  2,  where  a  =  1.674  and  e  =  0.1.  The  attracting  point  at  (a,/3)  = 
(0, 0)  is  removed  for  finite  €.  Trajectories  tend  to  a  <  0  where  no  collapse 
takes  place. 

Fig.  2  Results  of  numerical  simulations:  The  full  (open)  circles  correspond 
to  those  simulations  below  (above)  Pth-  The  solid  line  displays  the  analytic 
estimate  (28)  for  Pth  as  a  function  the  dispersion  parameter  7.  The  line, 

( — displays  the  estimate  using  linear  spreading,  and  the  dotted  line  using 
5  times  the  linear  estimate. 

Fig.  3  Focusing  geometries  considered  in  the  text:  (a)  The  incident  laser 
pulse  is  focused  from  air  through  a  lens  with  /  =  1.7  cm  directly  into  a 
water  sample,  where  it  comes  to  a  linear  focus  at  d  =  1.69  cm,  and  (b)  the 
incident  laser  pulse  is  focused  from  air  through  a  lens  with  /  =  1.7  cm,  back 
into  air  and  onto  a  thin  water  sample  at  the  linear  focus. 

Fig.  4  (a)  The  breakdown  intensity  Ibd  versus  pulse  duration,  and  (b)  the 
breakdown  input  power  Pbd  normalized  to  the  critical  power  for  collapse. 
Per  =  1  MW.  The  three  curves  correspond  to  initial  background  densities 
of  /jfc  =  0  (dotted  line),  pi,  =  10^°  cm”^  (solid  line),  and  pi,  =  10^^  cm“^ 
(dashed  line). 

Fig.  5  LIB  thresholds  with  SF:  (a)  Threshold  power  normalized  to  the 
critical  power  for  collapse,  and  (b)  the  threshold  intensity,  both  as  a  function 
of  pulse  duration  Tp.  In  (a)  the  solid  circles  give  the  threshold  for  a  peak 
density  of  10^®  cm“®,  the  triangles  are  the  threshold  powers  obtained  using 
the  criterion  of  a  peak  density  of  10^®  cm“®,  the  numbers  above  the  triangles 
are  the  %  absorption  of  energy  from  the  pulse  using  this  density,  and  the 
solid  curve  is  the  threshold  power  curve  obtained  without  SF,  as  in  Fig. 
2(b). 

Fig.  6  Peak  plasma  density  (in  units  of  10^®  cm“®)  and  percentage  of  energy 
absorption  as  functions  of  input  power  for  Tp  =  200  ps  (solid  line),  100  ps 
(dotted)  and  20  ps  (dashed),  with  pb  =  10^°  cm“®. 

Fig.  7  On-axis  maximum  field  maxt  |e(r  =  0,z,t)|  (top),  plasma  density 
maxt/>(r  =  0,z,t)  (middle)  and  total  pulse  energy  (bottom,  in  arbitratry 
units)  as  functions  of  propagation  distance  z  for  Tp  =  200  fs,  wq  =  120  pm 
and  P  =  l.OSFcr-  \ 
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Fig.  8  Surface  plots  of  |5(a:,0,z,t)|  (bottom)  and  p(x,0,z,t)  (top)  at  three 
propagation  distances:  (a)  z  =  1.687cm,  (b)  z  =  1.7cm  and  (c)  z  —  1.736cm. 

Fig.  9  Calculated  field  profiles  (solid  lines)  for  Eq  =  1.6  X  10®  V/m  and 
(a)  t=1.68  ps  before  the  field  enters  the  interface  at  zq  —  400  fim,  and  (b) 
t=2.75  ps  following  refiection  from  the  interface.  In  plots  (a)  and  (b)  we 
show  only  the  envelope  obtained  from  joining  the  field  peaks.  The  field 
strength  is  associated  with  the  left-hand  scale  and  is  normalized  to  unity. 
The  dashed  lines  in  (a)  and  (b)  show  the  corresponding  local  wavelength  over 
the  pulse,  and  are  associated  with  the  right-hand  scale.  Fig.  1  (c)  shows 
the  corresponding  pulse  spectra.  The  input  spectrum  is  shown  dashed  and 
is  associated  with  the  left-hand  scale. 

Fig.  10  Spatial  distribution  of  a)  the  full  field  in  the  vicinity  of  the  interface 
at  Zq  =  400  /xm,  and  b)  the  population  difference  n  at  times  t  =  1.61  ps 
(solid  line),  1.95  ps  (dashed  line),  2.28  ps  (single-dot-dashed  line),  and  2.62 
ps  (triple-dot-dashed  line),  where  the  front  propagates  to  the  right  with 
increasing  time,  and  c)  the  numerically  calculated  absorption  front  velocity 
normalized  to  c. 
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Abstract —  A  comprehensive  model  is  developed  for  focused 
pulse  propagation  in  water.  The  model  incorporates  self-focusing, 
group  velocity  dispersion,  and  laser-induced  breakdown  in  which 
an  electron  plasma  is  generated  via  cascade  and  muldphoton 
ionization  processes.  The  laser-induced  breakdown  is  studied  first 
without  considering  self-focusing  to  give  a  breakdown  threshold 
of  the  light  intensity,  which  compares  favorably  with  existing 
experimental  results.  The  simple  study  also  yields  the  threshold 
dependence  on  pulse  duration  and  input  spot  size,  thus  providing 
a  framework  to  view  the  results  of  numerical  simulations  of  the 
full  model.  The  simulations  establish  the  breakdown  threshold  in 
input  power  and  reveal  qualitatively  different  behavior  for  pico- 
and  femto-second  pulses.  For  longer  pulses,  the  cascade  process 
provides  the  breakdown  mechanism,  while  for  shorter  pulses 
the  cooperation  between  the  self-focusing  and  the  muldphoton 
plasma  generation  dominates  the  breakdown  threshold. 


1.  Introduction 

The  nonlinear  optical  properties  of  liquid  water  are 
of  considerable  interest  [\],  [2]  due  to  their  implications 
for  pulse  propagation  in  the  human  eye,  the  vitreous  humor 
being  predominantly  water.  Water  displays  a  nonlinear  Kerr 
effect  [2],  [3],  and  self-focusing  (SF)  can  occur  at  megawatt 
peak  input  powers.  In  addition  to  SF,  water  also  displays  laser- 
induced  breakdown  (LIB)  [4]-[6]  in  which  an  electron  plasma 
is  generated  due  to  the  high  intensity  of  focused  laser  pulses. 
The  electron  density  grows  explosively  and  serves  to  absorb 
and  scatter  the  remaining  pulse  energy,  thus  shielding  the 
area  beyond  the  focus  from  further  radiation,  i.e.,  the  retina 
[7]-[l0].  Indeed,  SF  leading  to  LIB  has  been  suggested  as  the 
explanation  for  anomolies  in  the  retinal  damage  data  for  visible 
femtosecond  pulses  [1 1]-I13].  Nonlinear  pulse  propagation  in 
water  is  therefore  of  considerable  interest  in  determining  laser 
safety  standards  in  the  ultrashort  pulse  regime. 
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In  this  paper,  we  report  numerical  results  on  linear  and 
nonlinear  propagation  in  liquid  water  using  200-ps  to  100-fs 
duration  pulses.  We  consider  a  geometry  in  which  a  Gaussian 
input  field  is  focused  through  a  lens  into  the  water  sample, 
thus  simulating  to  first  order  the  focusing  properties  of  the 
eye.  The  studies  incorporate  LIB,  SF,  multiphoton  absorption 
(MPA),  and  group  velocity  dispersion  (GVD).  The  model  is  an 
extended  nonlinear  Schrodinger  equation  (NLSE)  coupled  with 
an  equation  describing  the  plasma  generation.  The  NLSE  is 
widely  used  to  describe  light  propagation  and  has  been  studied 
extensively  [14].  An  interesting  mechanism  modeled  by  the 
NLSE  is  self-focusing  and  beam  collapse.  When  a  focused  or 
unfocused  propagating  laser  beam  undergoes  self-focusing,  the 
on-axis  intensity  gets  enhanced  considerably  if  the  power  of 
the  beam  is  near  but  below  a  critical  value.  When  the  power  is 
above  the  critical  value,  the  collapse,  namely  an  infinitely  high 
intensity,  occurs  at  finite  propagation  distance.  The  collapse 
is  clearly  unphysical  and  one  would  expect  other  nonlinear 
effects,  such  as  material  breakdown,  to  come  into  play  due  to 
the  high  intensity  produced.  Here,  we  extend  the  NLSE  to  in¬ 
clude  plasma  generation,  pulse-plasma  interaction,  and  MPA. 
The  plasma  generation  is  described  by  two  processes:  cascade 
(or  avalanche)  ionization  and  multiphoton  ionization.  Included 
also  is  the  GVD,  which  has  been  shown  to  lead  to  pulse- 
spliting  and  to  arrest  the  collapse  of  subpicosecond  pulses 
[i5]-[21].  Two  focusing  mechanisms  exist  in  our  studies: 
linear  focusing,  due  to  the  lens,  and  nonlinear  self-focusing. 
A  framework  for  these  studies  is  provided  by  first  studying 
LIB  without  SF,  and  this  model  compares  favorably  with 
existing  experimental  results.  We  then  establish  a  threshold 
for  LIB  by  numerical  simulations.  Our  results  show  that  SF  is 
a  major  player  in  understanding  LIB  in  water  for  pulses  of  the 
durations  considered  in  this  paper  and  of  spot-sizes  that  can 
be  easily  determined  from  a  criterion  given  in  this  paper. 

The  remainder  of  this  paper  is  organized  as  follows.  In 
Section  11,  we  describe  the  model  equations,  the  focusing 
geometries  and  material  parameters  pertinent  to  experiments, 
and  write ,  the  equations  into  dimensionless  form  with  ap¬ 
propriate  balings.  Although  we  will  present  our  results  in 
subsequent  sections  in  physical  units  to  facilitate  comparison 
with  experiments,  a  quick  inspection  of  the  dimensionless 
equations  leads  to  an  interesting  and  important  result,  i.e., 
the  relation  of  the  breakdown  threshold  intensity  to  the  pulse 
duration  for  the  case  where  multiphoton  ionization  is  the 
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The  threshold  at  which  self-focusing  initially  dominates  the  dynamics  of  short-pulse  propagation  in  normally 
dispersive  bulk  media,  causing  an  explosive  increase  in  peak  intensity,  is  estimated  analytically  and  verified 
numerically.  Intensity-dependent  propagation  effects  such  as  spectral  broadening  also  occur  explosively  at  this 
threshold. 


It  was  recently  proposed  by  Strickland  and  Corkum^ 
that  normal  group-velocity  dispersion  (GVD)  in¬ 
creases  the  peak  input  power  required  for  two- 
dimensional  (2D)  self-focusing  of  short  light  pulses. 
Their  proposal  is  intriguing  since  it  provides  a  lim¬ 
iting  mechanism  that  does  not  require  higher-order 
nonlinear  effects,  such  as  saturation  of  the  nonlinear 
refractive  index, ^  optical  breakdown,®  and  violation 
of  the  paraxial  approximation/ 

Because  cross  sections  of  a  light  pulse  differ  in 
power,  those  having  powers  that  exceed  the  critical 
power  self-focus  at  diferent  rates.  In  this  way,  self- 
focusing  tends  to  reduce  the  temporal  scale  length  of  a 
pulse,®'®  with  concomitant  increase  in  its  bandwidth, 
permitting  GVD  to  play  an  increasingly  important 
role  even  when  it  is  initially  a  small  effect.^ 

Normal  GVD  tends  to  ease  self-focusing  by  spread¬ 
ing  the  pulse  along  the  propagation  direction.  Re¬ 
cent  numerical  studies  of  the  nonlinear  Schrodinger 
equation  by  Chemev  and  Petrov®  and  Rothenberg® 
demonstrate  that  normal  GVD  increases  the  thresh¬ 
old  for  catastrophic  2D  self-focusing  through  a  tem¬ 
poral  pulse-splitting  process.  Analytic  research  has 
shown  that  the  critical  self-similar  collapse  singu¬ 
larity  associated  with  the  2D  nonlinear  Schrodinger 
equation  is  removed  by  normal  dispersion.^® 

The  simplest  and  most  ubiquitous  model  of  the 
spatiotemporal  propagation  of  a  light  pulse  is  ob¬ 
tained  as  a  uniform  as3rmptotic  expansion  of 
Maxwell’s  equations.  If  we  take  the  nonlinear  re¬ 
sponse  of  the  medium  to  be  instantaneous  and  cen- 
trosymmetric,  the  nonlinear  Schrodinger  equation 
for  the  envelope  of  a  light  pulse  is  obtained  as^  ® 


2ik 


-f- 


d(o  dt) 


+  vM  -  k 


d^k  d^A 
d(o^  dt^ 


+  2**^|APA  =  0. 

no 


(1) 


We  have  conducted  extensive  numerical  simulations 
of  Eq.  (1)  with  ^2  >  0,  d^kfdcj^  >  0,  and  initial  data 
of  the  form  A(x,y,2:  =  0,t)  =  exp[-(x^  +  y^)/2M;o^  - 
t^l2tp^]  to  explore  the  dynamics  of  pulse  propagation. 


Results  with  the  split-step  method^^  were  checked 
against  previous  numerical  simulations.®®’^^  The 
first  three  invariants  of  motion  and  the  magnitude 
of  higher-order  corrections  were  monitored  carefully. 
In  the  simulations  discussed  here,  nonparaxial  ef¬ 
fects  remained  small.  These  numerical  experiments 
show  that,  as  the  peak  input  power  is  increased,  a 
threshold  Pth  is  reached®  where  self-focusing  dom¬ 
inates  the  initial  dynamical  evolution  of  the  pulse. 
Near  Pth  the  pulse-splitting  process  coincides  with 
the  arrest  of  the  collapse  at  the  peak  of  a  pulse.  In 
this  Letter  it  is  shown  that  Pth  corresponds  to  the  in¬ 
put  power  at  which  the  length  scale  for  self-focusing 
equals  the  length  scale  for  nonlinear  dispersion.  In 
addition,  an  expression  for  Pth  is  derived  analytically 
and  verified  numerically. 

A  consequence  of  the  sharp  increase  in  peak  in¬ 
tensity  at  threshold  is  the  observation  of  intensity- 
dependent  processes.  An  important  example  is  the 
explosive  spectral  broadening  that  occurs  near  Pth- 
Through  four-wave  coupling  of  the  wave  trains  that 
make  up  the  pulse,  a  nonlinear  chirp  is  formed  along 
the  pulse.  As  a  result,  upshifbed  and  downshifted 
wave  packets  form  under  the  pulse.  The  two  wave 
packets  eventually  separate  in  the  presence  of  nor¬ 
mal  dispersion,  splitting  the  pulse.  The  spectrum  of 
the  resulting  double-peaked  pulse  is  broadband  and 
modulated.^*®  This  process  occurs  over  a  short  prop¬ 
agation  length  near  the  self-focusing  threshold,  re¬ 
sulting  in  an  explosive  increase  of  the  bandwidth. 
The  bandwidth  increases  more  than  would  be  ex¬ 
pected  by  self-phase  modulation  in  time.  The  gain  in 
the  four-wave  interaction  that  is  responsible  for  the 
broadening  is  favored  for  off-axis,  frequency-shifted 
modes.  Neither  pulse  splitting  nor  this  enhance¬ 
ment  in  the  spectral  broadening  is  captured  by  either 
stationary  self-focusing  or  one-dimensional  propaga¬ 
tion  alone. 

The  solutions  of  Eq.  (1)  are  functions  of  only  two  pa¬ 
rameters  when  Gaussian  initial  data  are  used.  We 
have  chosen  to  write  the  first  parameter  as  y  = 
which  measures  the  strength  of  the  dis¬ 
persion  relative  to  diffraction,  where  Lds  =  tp^jk” 
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The  SchrSdinger  equation  predicts  conical  emission  that  is  due  to  spatiotemporal  propagation  of  short 

pulses  in  normally  dispersive,  cubically  nonlinear  media.  This  effect  is  a  direct  consequence  of  a  four-wave 
interaction. 


Conical  emission  has  long  been  associated  with  spec¬ 
tral  superbroadening^  (SSB)  and  has  been  attributed 
to  four-photon  coupling.®  Recent  laboratory®  ex¬ 
periments  demonstrate  that  SSB  correlates  with 
the  threshold  for  self-focusing  (SF).  Both  these 
laboratory  experiments  and  several  numerical 
experiments*^  suggest  that  SF  is  inhibited  by 
normal  group-velocity  dispersion  (NGVD),  and  the 
threshold  required  for  the  observation  of  strong  SF 
effects  is  increased,®-®  having  a  simple  analytical 
form.®  Near  this  threshold  SSB,  conical  emission, 
and  pulse  splitting®  occur  explosively.  Here  we 
show  that  conical  emission  is  a  direct  consequence 
of  the  spatiotemporal  dynamics  of  short  laser  pulses 
described  by  the  nonlinear  SchrSdinger  equation 


2ik 


y  dz  d(o  dt ) 


+  Vx®A  -  k 


d^k  d®A 
dto^  df® 


-t-2ife®^|A|®A  =  0, 
no 


(1) 


where  A"  =  d^k/da)^^o  >  0  for  NGVD.  The  coni¬ 
cal  emission  predicted  by  Eq.  (1)  is  due  to  a  fun¬ 
damental  four-wave  interaction®  that  promotes  the 
transport  of  energy  to  a  band  of  modes  with  finite 
frequency  and  wave-number  shifts.  These  modes 
are  closely  related  to  those  of  the  modulational  insta¬ 
bility  of  the  condensate  or  the  uniform  plane-wave 
solution.^®  Pulse  splitting,  spectral  broadening,  and 
conical  emission  are  intimately  related,  and  each  is 
a  consequence  of  this  wave  interaction. 

For  short  pulses,  corrections  to  Eq.  (1)  appear  at 
the  next  order  in  the  uniform  asymptotic  expan¬ 
sion  and  accoimt  for  self-steepening,  stimulated 
Raman  scattering,  exponential  finite-time  nonlin¬ 
ear  response,^^  higher-order  dispersion,  nonparaxial 
effects,^^  and  vector  contributions.  Terms  in  the 
next  order  beyond  those  appearing  in  Eq.  (1)  intro¬ 
duce  asymmetries  in  time  and  in  frequency  space. 
Below  the  threshold  reported  in  Ref.  7  these  effects 
introduce  small  corrections  that  become  important 
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only  at  propagation  distances  larger  than  those  re¬ 
quired  for  pulse  splitting  and  conical  emission.  We 
neglect  higher-order  effects  and  focus  on  the  charac¬ 
teristics  properties  of  the  conical  emission  captured 
by  Eq.  (1). 

First  we  demonstrate  the  growth  of  frequency- 
shifted,  off-axis  modes  numerically  at  parameter 
values  below  the  threshold  reported  in  Ref.  7  and 
above  the  stationary  SF  threshold.  Pc-  For  the  pa¬ 
rameters  of  interest  here,  the  space  and  time  dy¬ 
namics  are  strongly  coupled,  and  standsu^d  anal3rtic 
technique  fail.  The  spUt-step  method^®  is  used 
to  integrate  Eq.  (1)  with  initial  data  of  the  form 
A{x,y,z  =  0,^)  =  exp[-{x^  +  y^)/wo^  -  {t/tp)^l  where 
Wo  is  the  beam  width  and  tp  is  the  pulse  length.  The 
parameters  were  chosen  so  that  Ldf  k''/tp^  —  0.1  and 
P/Pc  =  1.5,  where  Ldf  =  Wo^kl2  and  Pc  is  the  crit¬ 
ical  power  for  two-dimensional  SF.  In  Fig.  1  the 
contours  of  the  intensity,  |A(x,y  =  0,  f)P,  and  the  far- 
field  spectrum,  log[\A(kx,ky  —  0,  (y)|],  are  displayed 
for  z  =  0.0,  1.1,  and  1.5  times  the  diffraction  length, 
Ldf*  Figure  1(a)  shows  the  initial  pulse  as  it  enters 
the  medium.  Its  spectrum  is  locahzed,  appearing 
oval  in  (&> ,  k^).  In  Fig.  1(b)  the  pulse  has  propagated 
into  the  me^um,  reaching  a  stage  just  before  the 
splitting  event.  A  four-wave  interaction  has  trans¬ 
ported  energy  from  the  band  of  wave  trains,  (A,  w), 
to  sidebands  at  {kx  ±  k,,  ky  ±  Ky^w  ±  fl).  The  inter¬ 
action  broadens  the  spectrum  in  both  frequency  and 
wave  number,  but  it  is  most  efficient  along  cxuves 
for  which  where  =  (#Cx,Xy).  The  four- 

wave  interaction  that  produces  this  new  spectral 
featuie  dominates  the  spectral  broadening  process. 
In  Fig.  1(c)  the  pulse  has  split.  A  faint  indication  of 
modulations  in  the  frequency  spectrum  can  be  seen, 
which  at  higher  values  of  the  initial  power  are  read¬ 
ily  apparent.®  Analogous  spatiotemporal  four-wave 
interactions  destabilize  both  the  condensate^®  and 
the  cw  beam^^  solutions  of  Eq.  (1). 

A  consequence  of  the  growth  of  this  new  spectral 
feature  is  the  appearance  of  frequency-shifted  off-axis 
radiation  in  the  far-field  image  of  the  pulse,  which 
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.^umencal  integration  of  Maxwell’s  equations  for  propagation  of  a  femtosecond  pulse  in  a  medium 
with  linear  Lorentz  response  and  a  Kerr  nonlinearits  shows  shock  formation  on  the  underivinc  carrier 
‘he  envelope  shock.  The  earner  shock  is  characterized  by  the  appearance  of  a  strons 
third  harmonic  pu  se.  whereas  the  envelope  shock  appears  later  as  spectral  broadening  and  modulation 
or  tnc  lundamental  and  higher  harmonic  spectral  features. 


P.ACS  numbers:  42.65. Ky 


Advances  in  laser  technology  in  the  past  decade  have 
made  possible  the  production  of  pulses  which  contain  a  few 
optical  cycles  1 1  ].  Although  such  ultrashort  pulses  contain 
small  amounts  of  optical  energy,  enormous  intensities 
exceeding  1  TW/cm‘  can  arise  and  the  accompanying 
intensity-dependent  corrections  to  the  index  of  refraction 
are  such  that  one  can  expect  novel  nonlinear  phenomena 
such  as  shock  formation  over  very  short  propagation 
lengths.  It  has  been  suggested  in  \2]  that  an  envelope 
shock  was  observed  experimentally  [3]  and  that  this  can  be 
understood  using  standard  envelope  approximations  (2.4] 
to  Maxwell's  equations.  Laser-induced  breakdown  (LIB) 
cannot  be  ruled  out  at  such  very  high  peak  intensities, 
but  there  is  evidence  to  show  that  for  such  short  and 
hence  low  energy  pulses  the  cascade-avalanche  path  is 
unlikely  and  multiphoton  processes  are  more  likely  to  lead 
to  breakdown  [5].  Indeed,  very^  recent  expieriments  in 
water  using  100  fs  pulses  indicate  that  local  peak  field 
intensities  can  exceed  10'^  W/cm’  in  the  focal  region 
with  incident  pulse  absorption  being  less  than  5%  [6]. 
Moreover,  breakdown  becomes  a  sensitive  function  of 
optical  wavelength.  Shock  formation  on  the  carrier  wave 
is  expected  therefore  to  compete  with  other  physics  during 
the  critical  collapse  of  femtosecond  duration  optical  pulses 
in  optically  transparent  media  where  the  local  intensity  at 
the  critical  collapse  distance  can  become  very  large. 

The  above  breakdown  scenario  is  extremely  compli¬ 
cated,  so  we  confine  our  attention  here  to  plane  wave 
propagation  for  simplicity  and  show  that  an  optical  carrier 
shock  can  arise  in  a  medium  with  an  instantaneous  Kerr 
nonlinearity.  Dispersion  plays  an  important  role  in  shock 
regularization  (smoothing)  and  influences  the  signature  of 
the  c^er  shock.  As  dispersion  is  typically  strong  for  such 
short  optical  pulses  [the  dispersion  length  scales  as 
where  k"  is  the  leading  order  contribution  to  the  group  ve¬ 
locity  dispersion  (GVD)  and  Tp  is  the  characteristic  pulse 
length],  phase  mismatch  leads  to  the  separation  in  time  of  a 
strong  third  harmonic  optical  pulse  moving  with  a  different 
group  velocity  from  the  fundamental.  For  very  weak  dis¬ 
persion,  a  component  of  the  third  harmonic  pulse  moves 
with  the  fundamental  and,  in  the  dispersionless  case,  all 
higher  harmonics  of  the  fundamental  are  phase  matched 


and  see  explosive  growth.  .As  the  phenomena  we  are  con¬ 
cerned  with  occurs  on  the  scale  of  the  carrier  w  avelencth. 
no  envelope  approximations  are  valid  and  one  must  resort 
to  a  numerical  integration  of  .Maxwell's  equations.  .Nu¬ 
merical  schemes  for  the  integration  of  Maxwell's  equations 
have  been  refined  in  the  past  few  years  [7]  to  allow  for  an 
efficient  integration  of  media  with  memory  in  both  the  lin¬ 
ear  and  nonlinear  polarizations. 

We  restrict  our  attention  to  nonmagnetic  dielectric  me¬ 
dia  with  no  free  charges,  in  which  case  we  have  for 
Maxwell's  equations 

hz  bt  ft-i)  bz  ' 

where  all  quantities  above  and  in  the  following  are  in 
MKS  units.  The  medium  is  modeled  by  a  single  Lorentz 
oscillator  plus  an  instantaneous  Kerr  nonlinearity 

Dj(c.r)  =  eo|ex£,(c.r)  -  ^  dt' ~  t')E^{z..t') 

with  the  linear  susceptibility  given  by  xU)  =  -  X 

sin(V<uj5  -  5*/4r)/-y<u5  ixico)  =  - 

iSu>  —  tu')],  (Op  =  (e^  —  £i)£t>o.  Cj  and  e*  are  the 
Static  and  infinite  relative  permittivities,  respectively,  and 
.  ioo  the  resonance  frequency  of  the  Lorentz  oscillators. 
Maxwell  s  equations  are  solved  by  either  a  second  order 
in  time,  second  order  in  space  [(2,2)]  finite  difference  time 
domain  method  [8]  or  a  second  order  in  time,  fourth  order 
in  space  [(2,4)]  scheme  [9].  The  numerical  dispersion 
inherent  in  these  methods  has  recently  received  a  good 
deal  of  attention  [10],  and  we  have  chosen  our  spatial 
discretizations  accordingly. 

The  convolution  integral  / dt'  xU  ~  t')Exit') 

in  Eq.  (2)  is  most  efficiently  solved  by  replacing  it  with 
the  equivalent  second  order  ordinary  differential  equation 
[7,11]  ^ 


a>o 


- L. 


uiq  dt 


+  Pl  =  — 2  €oE:,  , 

<tf0 


(4) 


- ^ 

which  is  solved  by  second  order  central  differencing 
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We  present  numerical  studies  of  nonlinear  propagation  for  picosecond  pulses  focused  in  water.  Depending 
on  the  pulse  duration  and  focusing  conditions,  for  some  input  powers  self-focusing  may  precede  laser-induced 
breakdown  and  \'ice  versa.  We  derive  a  criterion  that  predicts  the  relative  roles  of  laser-induced  breakdown 
and  self-focusing,  t  1995  Optical  Society  of  America 


The  nonlinear-optical  properties  of  liquid  water  are  of 
considerable  interest^-^  owing  to  their  implications  for 
propagation  in  the  eye,  the  vitreous  humor  being  pre¬ 
dominantly  water.  Water  displays  self-focusing  (SF) 
at  megawatt  input  powers^*^  and  has  been  suggested 
as  the  source  of  retinal  damage  at  anomolously  low 
energies  for  femtosecond  pulses.^  Water  also  displays 
laser-induced  breakdown  in  which  an  elec¬ 

tron  plasma  grows  exponentially  by  avalanche  ioniza¬ 
tion.'*  Once  generated,  this  electron  plasma  serves  to 
absorb  and  scatter  the  remaining  pulse  energy,  thus 
effectively  shielding  the  area  beyond  the  plasma,  i.e., 
the  retina. The  interplay  between  SF  and  LIB  is 
therefore  of  importance  in  assessing  laser  eye  damage. 

In  this  Letter  we  report  numerical  simulations  of 
SF  and  LIB  in  water  using  picosecond  pulses.  Our 
numerical  simulations  revealed  cases  of  peak  input 
powers  for  which  SF  preceded  LIB  and  vice  versa,  de¬ 
pending  on  the  pulse  duration  and  focusing  conditions 
(see  also  Refs.  2,  12,  and  13).  Here  we  develop  a  criti¬ 
cal  power  Plib  for  LIB  that  can  be  compared  with  the 
critical  power  Per  for  SF:  For  Plib  >  Pery  SF  precedes 
LIB;  for  Plib  <  Per,  LIB  precedes  SF.  Numerical 
simulations  are  presented  to  show  that  this  criterion 
predicts  the  correct  dependence  on  focal  length  and 
pulse  duration. 

Our  theoretical  model  for  LIB  is  that  previously 
described  by  Feit  and  Fleck^^  appropriately  extended  to 
include  nonlinear  SF.^^  The  equation  for  the  electric- 
field  envelope  P (r,  z,  t)  in  a  reference  frame  moving  at 
the  group  velocity  is  then 

^  -  yd  +  iwT)nT  +  ikon2\TfT. 

(1) 

along  with  the  following  Drude  model  describing 
avalanche  ionization  of  the  electron  density  n  (Ref.  6): 

-  On  2)  ^nirp- an*.  (2) 

Here  is  the  intensity,  k  =  =  nfcoi/c,  112  is  the 

nonlinear  coefficient,  <r  is  the  cross  section  for  electron- 

ni /f o o  1  o c:o .  AA 


neutral  inverse  bremsstrahlung,  t  is  the  electron 
collision  relaxation  time,  Eg  is  the  ionization  energy, 
and  a  is  the  rate  of  radiative  recombination.  Here  we 
consider  the  following  parameter  values  appropriate 
to  the  visible  wavelength  Aq  =  580  nm:  r  =  10"^®  s,^^ 
rib  =  1-33,  n2  =  2  X  10" cm^/W,^  Eg  =  6.5  eV,'' 
a  -  10^®  cm^  s'^*'  and  <7  =  1.4  x  10"^'  cm^.^^ 
The  critical  power  for  SF  is  Per  =  0.15A2/n6n2,  which 
yields  Per  =  2  for  the  above  parameters.^ 

We  are  interested  in  the  solutions  of  Eqs.  (1)  and  (2) 
for  an  initial  Gaussian  beam  focused  onto  the  water 
sample  through  a  lens  of  focal  length  f. 


w^here  Pjn  is  the  peak  power,  wq  is  the  spot  size,  and  tp 
is  the  pulse  length.  In  the  absence  of  LIB  the  behavior 
of  the  field  is  dictated  by  the  power  ratio  Pin /Per  since 
each  time  slice  of  the  pulse  acts  independently.^®  For 
Pin /Per  >  1  the  central  time  slice  of  the  pulse  undergoes 
collapse  to  a  singularity  at  a  nonlinear  focal  distance 


—  =  _  H-  - 1 - ,  (4) 

fNL  f  ZNL(P)n) 

where  ZNL(Pin)  is  the  nonlinear  collapse  distance  with¬ 
out  the  lens.^®  For  Pm/Per  =  1  we  have  znl 
that  the  field  collapses  at  the  linear  focus. 

Our  numerical  simulations  based  on  Eqs.  (1)  and  (2) 
show  that  LIB  prevents  catastrophic  SF  to  a  singular¬ 
ity  by  absorption  and  defocusing  of  the  incident  field 
that  result,  from  the  electron  plasma.  However,  SF 
can  proceed  a  considerable  distance  into  the  collapse 
before  it  is  arrested  by  LIB.  To  obtain  a  criterion  for 
the  relative  roles  of  SF  and  LIB  we  neglect  the  effects 
of  SF.  Then  the  input  Gaussian  field  wdll  come  to  a 
focus  at  a  longitudinal  position  d  =  /*/(!  +  f^/zo)^ 
where  z©  “  fl'Wo^no/Ao  is  the  Rayleigh  range  of  the  in' 
put  beam.  The  size  of  the  laser  beam  Wf  at  the  focus 
is  (Wf/wo)^  ^  with  Zf  th^ 

Rayleigh  range  of  the  focused  beam.  To  obtain  a  cri* 
terion  for  LIB  we  use  the  following  model:  The  laser 
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We  calculate  the  spatiotemporal  evolution  of  intense,  femtosecond  pulses  that  are  incident  upon  a  saturable 
absorbing  interface  in  the  regime  of  self-reflection,  using  the  finite-difference  time-domain  computational 
method.  The  pulses  induce  a  curved,  moving  absorption  front  in  the  nonlinear  medium  that  acts  as  a  transient 
focusing  mirror  for  the  reflected  pulse  energy.  ©  1997  Optical  Society  of  America 


The  propagation  of  an  intense  optical  field  that  is  inci¬ 
dent  from  air  upon  a  saturable  absorber  can  give  rise 
to  a  self-reflected  wave  if  the  absorption  is  sufficiently 
large.^  The  reflected  wave  arises  at  a  wavelength- 
scale  spatial  transition  between  the  saturated  and  the 
unsaturated  regions  of  the  absorber,^  akin  to  the  skin 
effect  at  a  linear  boundary,  and  is  self-reflected  in 
the  sense  that  the  transition  region  is  induced  by  the 
incident  field  itself.  The  self-reflected  wave  first  pre¬ 
dicted  for  a  cw  incident  plane-wave  field'  was  subse¬ 
quently  extended  to  consideration  of  a  cw  incident  beam 
with  a  Gaussian  transverse  profile.®  In  the  latter  case 
the  emergence  of  a  halo  was  predicted  in  the  far  field  at 
high  intensities.  Recently,  the  finite-difference  time- 
domain  (FDTD)  method^’^  was  used  to  study  the  dy¬ 
namics  of  self-reflection  under  pulsed  excitation  and  to 
predict  an  intensity-dependent  Doppler  shift  in  the  re¬ 
flected  pulse.® 

In  this  Letter  we  include  transverse  variations  to 
study  the  spatiotemporal  dynamics  of  self-reflection. 
The  FDTD  method  discretizes  the  differential  form 
of  Maxwell’s  equations  directly  and  allows  one  to  de¬ 
termine  the  evolution  of  an  optical  field  in  a  nonlin¬ 
ear  medium,  subject  to  the  given  constitutive  relations 
between  the  electric  field  and  polarization  (see,  e.g., 
Refs.  7  and  8)  and  without  recourse  to  the  slowly  vary¬ 
ing  envelope  approximation.  We  examine  pulse  shap¬ 
ing  in  the  near  field,  close  to  the  absorber  boundary, 
and  predict  a  new  nonlinear-focusing  effect,  which  we 
attribute  to  the  formation  of  a  transient  focusing  mir¬ 
ror  in  the  absorber.  As  the  incident  pulse  impinges 
upon  and  strongly  saturates  the  absorber,  it  excites 
a  moving  reflection  front®  that  is  shaped  according 
to  the  transverse  profile  of  the  incident  pulse.  If  the 
incident  pulse’s  transverse  profile  is  bell  shaped,  so 
too  is  the  resulting  mirror,  and  the  reflected  pulse 
is  focused  according  to  the  waist  and  intensity  of  the 
incident  pulse.  In  addition  to  being  reshaped,  the 
reflected  pulse  is  spectrally  broadened  and  red  shifted 
owing  to  the  Doppler  effect  at  the  moving  mirror.® 

We  consider  the  time-dependent  propagation  of  a 
two-dimensional  TE-polarized  pulse,  in  which  the  elec¬ 
tric  field  is  polarized  along  the  y  axis  and  is  also 
assumed  to  be  uniform  along  that  axis,  E(r,0  — 
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yEy{x,  z,  t).  Then  the  Maxwell  equations  for  the  elec¬ 
tric  and  magnetic  field  quantities  Eyj  Bx,  and  Bz  are 

1  (^Bx  dBz\ 

dt  dz'  dt  dx'  dt  fio\  dz  ax  / 

(1) 

where  the  z  axis  is  the  propagation  direction  and  the  x 
axis  is  the  transverse  direction.  The  nonlinear-optical 
response  of  the  saturable  absorber  is  included  by  use 
of  a  two-level  model  through  the  constitutive  relation 
Dy  =  €oEy  +  Py,  whcrc  the  macroscopic  polarization 
Py  ==  Np{p2i  +  c.c.)  is  determined  by  the  Bloch  equa¬ 
tions® 

+  (  72  +  io)2l)p2\  -  i 

^  +  7i(n  -  1)  =  2i  ^^(p2i  ““  P21*)*  (2) 

Here  N  is  the  density,  p2i  is  the  off-diagonal  density 
matrix  element,  n  =  (pn  -  P22)  is  the  population 
difference,  <1^21  is  the  transition  frequency,  p  is  the 
dipole  moment,  and  yi  =  1/ti  and  72  =  1/^2  are 
the  population  and  polarization  damping  constants, 
respectively. 

The  conditions  for  self-reflection  of  a  continuous 
plane-wave  incident  field  require  that  the  normalized 
parameters  t//  =  p^N/60^72  and  F  =  p£o/^(7i72)^^ 
be  greater  than  unity,  with  £0  the  peak  input  field. 
Physically,  this  requires  that  the  linear  absorption  be 
large,  on  a  wavelength  scale,  and  that  the  incident 
field  be  strong  enough  to  saturate  the  absorption.^  In 
the  case  of  ultrashort  pulses  we  also  require  that 
the  incident  pulse  duration  tp  be  greater  than  the 
polarization  decay  time  but  less  than  the  population 
decay  time,Ui  >  tp  >  ^2•^  To  meet  these  conditions 
for  the  sub-lOO-fs  pulses  to  which  we  were  restricted 
by  computational  resources,  we  adopted  the  following 
medium  parameters  for  illustrative  purposes:  ti  = 
0.5  ns,  t2  =  10  fs,  (o  =  ct>2i  =  2  X  10'^  rad  s"'  (A  = 
942  nm),  p  =  4  x  lO'^^  C,  and  No  =  2  x  lO^®  cm“3. 
For  these  values  a  normalized  field  strength  of  F  = 
1  corresponds  to  an  electric  field  strength  of  £0  = 
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We  introduce  the  dynamic  nonlinear  optical  skin  effect  in  which  a  pulse  incident  on  a  saturable 
absorbing  interface  is  self-reflected  from  a  moving  absorption  front.  The  motion  of  the  front  causes 
the  self-reflected  wave  to  be  redshifted  by  the  Doppler  effect,  which  in  turn  serves  as  an  experimentally 
observable  signature  for  the  front  propagation.  [50031-9007(96)00202-5] 


PACS  numbers:  42.65 .-k,  41.20,Jb,  42.50.Gy 

In  the  linear  optical  skin  effect  a  pulse  incident  from 
air  is  reflected  from  a  highly  absorbing  interface  after 
penetrating  only  a  fraction  of  a  wavelength  into  the  ab¬ 
sorbing  medium,  this  distance  being  the  skin  depth  [1,2]. 
The  skin  effect  is  therefore  of  fundamental  importance  in 
understanding  the  electrodynamics  of  pulse  propagation 
at  condensed  matter  interfaces,  such  as  metals  for  field 
frequencies  below  the  plasma  frequency  [2],  and  semicon¬ 
ductors  with  highly  absorbing  excitonic  features  [3],  In 
addition,  it  belongs  to  an  important  class  of  optical  prob¬ 
lems  for  which  the  notion  of  an  electromagnetic  field  en¬ 
velope  varying  slowly  on  the  scale  of  a  wavelength  simply 
does  not  apply.  The  skin  effect  cannot  be  understood  on 
the  basis  of  envelope  equations  but  is  rather  a  consequence 
of  Maxwell’s  equations  for  the  interface. 

In  this  Letter  we  introduce  the  dynamic  nonlinear 
optical  skin  effect  for  pulses  and  elucidate  the  underlying 
physics.  In  the  nonlinear  skin  effect  a  high  intensity 
pulse  is  incident  upon  a  nonlinear  absorbing  interface. 
Broadly  speaking,  saturation  of  the  absorption  allows  the 
incident  field  to  penetrate  beyond  the  linear  skin  depth 
into  the  medium,  and  this  causes  an  absorption  front  to 
propagate  into  the  medium  which  separates  the  regions  of 
low  (saturated)  and  high  (unsaturated)  absorption.  The 
front  is  excited  by  the  incident  pulse  which  is  in  turn 
reflected  from  the  sharp  absorption  front,  yielding  a  self- 
reflected  pulse  [4].  Thus  the  absorption  front  acts  as  a 
moving  mirror  from  which  the  pulse  is  self-reflected,  and 
the  pulse  suffers  a  redshift  due  to  the  Doppler  effect  [5]. 

Continuous  wave  (cw)  self-reflection  from  stationary  ab¬ 
sorption  fronts  for  plane  wave  [4]  and  transverse  Gaussian 
[7]  fields  incident  at  sharp  and  smooth  [6]  interfaces  has 
been  studied  theoretically  but  not  experimentally  verified 
so  far.  In  part,  this  is  due  to  the  extremely  high  absorption 
and  strong  saturation  required  for  its  manifestation,  but  the 
difficulty  of  obtaining  good  experimental  signatures  should 
not  be  overlooked.  Here  we  explore  the  transient  regime 
using  the  two-level  Maxwell-Bloch  equations.  In  partic¬ 
ular,  we  show  that  moving  fronts  are  excited  by  the  inci¬ 
dent  pulse  [8,9]  and  that  the  self-reflected  pulse  bears  clear 


spectral  signatures  due  to  the  Doppler  effect,  which  should 
be  observable  experimentally. 

We  consider  the  time-dependent  propagation  of  a 
linearly  polarized  plane  electromagnetic  wave  incident  on 
a  nonlinear  medium  composed  of  two-level  systems.  For 
propagation  along  the  z  axis,  and  taking  the  electric  field 
polarized  along  the  x  axis.  Maxwell’s  curl  equations  take 
the  form  [1,2] 

dBy  _  ^  _  ClHy 

dt  hz  8t  hz 

where  By  =  fioHy.  The  specification  of  the  problem  is 
completed  with  the  constitutive  relation  Dy  =  6()£v  + 
Py,  where  Py  is  the  optical  polarization.  To  elucidate  the 
basic  physics  we  employ  a  two-level  model  to  describe  the 
optical  response  with  lower  electronic  state  |1)  and  upper 
state  |2).  The  Bloch  equations  are  then  (see,  for  example. 
Ref.  [10]) 

^P2\  /  .X  ^ 

— - (72  +  i<^2\)P2\  —  t  ^ 

ot  Ft 

^  +  7i(«  -  1)  =  2/ ^  (P2I  - />2i).  (2) 

where  p2i  is  the  off-diagonal  density  matrix  element,  n  = 
Pii  P22  is  the  population  difference  between  the  lower 
and  upper  states,  <02\  is  the  transition  frequency,  p  is  the 
dipole  moment  in  the  field  direction,  and  y\  and  72  are 
phenomenological  damping  constants  for  the  population 
and  polarization,  respectively.  The  polarization  due  to 
the  atoms  is  then  given  by  P,v  =  V(z)p(p2i  +  c.c.),  with 
N{z)  the  density  of  two-level  systems  which  varies  along 
z  in  general. 

Equations  (1)  and  (2)  are  solved  using  a  standard 
discretization  scheme  described  by  Yee  [11]  and  the 
Bloch  equations  integrated  in  time  using  a  fourth-order 
Runge-Kutta  method.  The  system  of  equations  was 
solved  with  the  initial  condition  on  the  field 

Ex{zJ  =  0)  =  £ocos[27rtu(z  - 

(3) 
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REFLECTION  OF  LOCALIZED  BEAMS  FROM  A  NONLINEAR 

ABSORBING  INTERFACE  * 

J.  A.  POWELL^  E.  M.  WRIGHT^  AND  J.  V.  MOLONEY§ 

Abstract.  Beam  propagation  and  reflection  are  studied  at  the  interface  between  a  Kerr  medium 
and  a  saturable  absorbing  nonlinear  medium.  For  spatial  solitons  propagating  nearly  tangent  to  the 
interface,  slowly  varying  envelope  theory  describes  beam  behavior.  Phenomena  are  studied  using 
numerical  beam  propagation  and  dynamical  systems  arguments.  Projecting  the  dynamics  onto  the 
sohton  modes  of  the  nonabsorbing  medium,  an  effective-particle  theory  is  developed  for  the  reflected 
beain.  The  results  of  the  effective-particle  theory  and  the  numerical  analysis  are  compared,  and  a 
heuristic  criterion  for  predicting  reflection  and  absorption  of  beams  is  derived. 

Key  words,  beam  reflection,  nonlinear  optics,  effective-particle  theory 
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1.  Introduction.  The  reflection  properties  of  nonlinear  dielectric  interfaces  have 
attracted  considerable  interest  over  the  last  years  due  to  potential  applications  in  all- 
optical  switching  [1]-[13].  The  basic  idea  is  to  understand  what  adjustable  parameters 
control  the  reflection  of  beams  at  interfaces  between  optical  materials.  The  angle  of 
beam  reflection  may  be  controlled  simply  by  varying  the  power  or  phase  of  an  input 
beam,  without  changing  the  incident  angle  of  the  input  beam.  This  could  result  in 
switching  devices  that  are  all-optical,  that  is,  that  do  not  depend  on  conversion  of 
optical  to  electric  signals,  and  subsequent  slowing  of  response  and  loss  of  information, 
for  switching.  In  particular,  using  plane  wave  theory,  Kaplan  first  predicted  bistable 
reflection,  or  the  possibility  of  multiple  beam  reflections  dependent  on  external  con¬ 
trol  parameters,  from  a  linear-nonlinear  interface  for  incidence  angles  close  to  the 
criti(^  angle  for  total  internal  reflection  [1].  Subsequent  numerical  simulations  using 
localized  beams  of  finite  transverse  extent,  typically  Gaussians,  failed  to  show  bistable 
reflection  but  exposed  several  exciting  new  phenomena  such  as  the  nonlinear  Goos- 
Hanchen  effect  [9],  and  the  transmission  of  self- focused  channels  (or  solitons)  through 
the  dielectric  interface  [2],  [7]. 

The  dynamics  of  localized  beams  incident  on  nonlinear  dielectric  interfaces  was 
put  on  a  firm  theoretical  foundation  by  the  introduction  of  the  effective-particle  theory 
by  Aceves,  Moloney,  and  Newell  [10].  In  this  theory  the  incident  beam,  usually  a  spatial 
soliton,  is  treated  as  a  particle  whose  motion  is  dictated  by  an  effective  potential  that  is 
derived  from  the  linear  and  nonlinear  properties  of  the  interface.  Thus,  for  example, 
the  nonlinear  Goos— Hanchen  shift  can  be  viewed  as  critical,  slowing  down  as  the 
particle  approaches  a  turning  point  of  the  effective  potential. 

In  this  paper  we  consider  the  reflection  of  localized  beams  from  nonlinear  ab¬ 
sorbing  interfaces.  Roso-Franco  previously  showed  that  for  plane  waves  at  normal 
iiicidence  the  nonlinear  absorbing  interface  can  display  many  novel  features  including 
bistability  and  the  generation  of  backward-propagating  waves  [14].  Here  we  investi- 
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